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SUMMARY 


A numerical simulation of bubble-type vortex breakdown using a unique 
discrete form of the full three-dimensional, unsteady incompressible Navier- 
Stokes equations was performed. The Navier-Stokes equations were written in 
a vorticity-velocity form and the physical problem was not restricted to 
axisymmetric flow. Based on the results of a previous study, the problem 
was parameterized in terms of a Rossby number-Reynolds number basis. Utili- 
zation of this parameter duo was shown to dictate the form of the free-field 
boundary condition specification and allowed control of axial breakdown 
location within the computational domain. The structure of the breakdown 
bubble was studied through time evolution plots of planar projected velocity 
vectors as well as through plots of particle traces and vortex lines. These 
results compared favorably with previous experimental studies. In addition, 
profiles of all three velocity components are presented at various axial 
stations and a Fourier analysis was performed to identify the dominant 
circumferential modes. The dynamics of the breakdown process were studied 
through plots of axial variation of rate of change of integrated total 
energy and rate of change of integrated enstrophy, as well as through 
contour plots of velocity, vorticity and pressure. 
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NOMENCLATURE 


an adjustable constant associated with the vortex core 
diameter 

n x m row vector where m is the number of unknowns in a 
cell 

3x3 exponential transforma tion matrix 

a vector function, the cross product of which is 
divegrence free 

. 3u. 3u- 

defined as j ( 3 - * 
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scalar quantity in Kaczmarz relaxation scheme 
boundary condition for the velocity vector 
defined as (Ax/2, Ay/2, Az/21 respectively 

computational domain lengths in the x, y and z directions, 
respectively 

a constant proportional to circulation 
pressure variable 
radial coordinate 

radial reference length defined as ( 2 v/a) ^ 2 

radius where the swirl velocity is a maximum 

residual quantity in Kaczmarz relaxation scheme 

★ ★ 

U r 

Reynolds number; defined as — - — 
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time variable 

velocity components in (x,y,z) system of coordinates 
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axial velocity at a radius equal to r* 
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axial velocity excess (or deficit) at the vortex 
centerl ine 

swirl velocity 

Cartesian system of coordinates: x-axial, y-transverse, 
z-spanwise 

nondimensional coordinates for use in grid generation; 
0<(x,y,z)<l 

transverse and spanwise coordinates corresponding to the 
location of the vortex centerline 

defined as (Ax/Az,Ax/Ay,Ay/Az) , respectively 

defined as U /U 
o • 

central difference spatial operators 
central difference time operator 
discretization intervals 

stretching parameter for grid in y and z directions 

transformed vorticity in (x,y,z) system of coordinates, 
respectively 

defined as (-j^, , respectively 

defined as (At/ Ax, At/ Ay, At/ Az), respectively 

spatial averaging operators 

time average operator 

kinematic viscosity 

den si ty 

stretching parameter for grid in x direction 

defined as t j, • t. 

n+ V 2 n 

defined as t n+1 - t n+ i /2 

variables related to the gradient of vorticity 

a scalar quantity used in the Helmholtz projection 

vorticity components in (x,y,z) system of coordinates 
respectively 

rate of rotation taken in the limit as r+o 
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INTRODUCTION 


Vortices and vortical motions have played an important role in the 
development of theoretical fluid mechanics. For instance, the Helmholtz 
Theorems of Vorticity and generalizations by Kelvin, Crocco and others 
established flow properties involving the kinematics of vortex lines and 
the dynamics of vorticity. The theory of lifting surfaces, developed by 
Prandtl, Kutta and Joukowski, is based on the concept of a bound 
vortex. Recently, the recognition of large scale coherent vortical 
structures in turbulent flows has resulted in renewed interest in the 
study of vortices. 

Although no general definition of a vortex exists, it can be 
thought of as a collection of fluid particles rotating around a common 
axis. Mathematically, vorticity is defined as the curl of the velocity 
vector and is equivalent to twice the angular velocity of a fluid 
particle. In addition, it is not necessary for a vorticity field to 
represent a vortex, an example being a parallel shear flow. 

The most common development of a vortex occurs when a boundary 
layer separates from a surface and rolls up into a wake vortex. Tip 
vortices fall into this class. Tip vortices are generated when a fluid 
flows against a finite plate or sharp edged body at a nonzero angle of 
attack. These vortices are often highly stable structures and are 
characterized by a strong axial flow. Other examples of vortices with 
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an axial velocity component include tornados and waterspouts, intake 
vortices and swirling flows in pipes and tubes. 

The presence of tip vortices in the wake of large aircraft 
constitute a hazard in areas of dense air traffic. These vortices can 
cause severe rolling of smaller aircraft that enter them. They are 
dissipated either by viscous dissipation, a sinusoidal type instability 
or, infrequently, by a core bursting mechanism. 

Leading edge vortices shed from a delta wing induce a velocity 
field that results in increased lift and stability of the wing [1]. 
However, under certain conditions related to the angle of attack of the 
wing, these vortices can undergo a sudden and drastic change in 
structure known as vortex breakdown. This breakdown can alter severely 
the aerodynamic characteristics of the wing. 

Swirling flows have been used to stabilize high intensity 
combustion processes [2]. Here a recirculation zone acts to stabilize 
combustion by recirculating hot gases to the root of the flame. In 
addition, combustion lengths are reduced due to the high levels of 
entrainment induced by the recirculation zone. 

The ability to control these vortical structures is an important 
and active area of research. For example, it is desirable to delay the 
breakdown process over a delta wing and accelerate the process in 
regards to trailing wing tip vortices. In combustion applications, the 
internal structure of the recirculation (breakdown) region is of 
critical importance. Unfortunately, a comprehensive theory to describe 
the breakdown process and the parameters affecting it is lacking 


presently, although several have been proposed. A review of these 

theories and supporting numerical and experimental work follows. 

Yortex breakdown was first observed experimentally by Peck ham and 
Atkinson [3], They observed that vortices shed from a delta wing at 
high angles of attack appeared to "bell out" and dissipate several core 
diameters downstream from the trailing edge of a wing. Since then, 
vortex breakdown has been observed in swirling flows in straight pipes, 
nozzles, diffusers and combustion chambers, [2,4] and tornados [5]. 
Seven types of breakdown have been identified experimentally, [6] 
ranging from a mild "spiral" type to a strong "bubble" type breakdown. 
Observations in the early 1960's spurred considerable effort to develop 
a theoretical explanation of the vortex breakdown phenomena. Three 

schools of thought can be identified, all of which may be divided into 
three separate groups: 1) The concept of a critical state [7,8,9]; 
2) Analogy to boundary layer separation [10,11] and, 3) Hydrodynamic 

instability [12,13,14]. 

Squire [7] appears to be the first to have performed a theoretical 
analysis of vortex breakdown. He suggested that if standing waves were 
able to exist on a vortex core then small disturbances, present 

downstream, could propogate upstream and cause breakdown. This is 
analogous to the earlier work of Taylor [15] on the stability of 
circular Couette flow. There, a linear stability analysis was performed 
to ascertain the ability of the base flow to support axisymmetric 

standing wave disturbances. Two of the cases studied by Squire assumed 
that the vortex flow was inviscid and axisymmetric. The assumed form of 
the upstream velocity distribution resulted in a linear disturbance 
equation which he then solved to determine a condition under which an 
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inviscid, a xi symmetric, steady perturbation to the flow could exist. 
This condition, which was necessary for the existence of a standing 
wave, was taken to mark the transition between subcritical and 
supercri tical states. . In Squire's first case, the axial velocity, U, 
was considered to be a constant. The dimensionless swirl velocity V was 
taken piecewise as that of solid body rotation inside a core of unit 
radius (V=V Q r) and connected with a potential vortex outside (V=V Q /r). 
A constant, V 0 , was used to control the swirl. He found that for 
standing waves to exist a swirl parameter, "k", which was the ratio of 
the maximum swirl speed to the axial speed (Y max /U), had to be greater 
than, or equal to 1.20. When k=1.20 the wave is infinitely long, but it 
has a finite wavelength for k>1.20. 

In the second case U was also taken to be a constant, but the swirl 

-r 2 

velocity was assumed to be V 3 (V Q /r) (1 — e ) with V 0 a nondimensional 

parameter. Again Squire found that there was a condition on the swirl 
parameter k for the existence of a standing wave. This condition was 
k 33 V [nax />1.0 where we note that ^max® 0 * 688 V 0 . 

Benjamin [8] examined this phenomena from a different point of 
view. He considered vortex breakdown to be a finite transition between 
two dynamically conjugate states of flow, similar to the occurrence of a 
hydraulic jump in open channel flow. These are a subcritical flow, 
which is defined as the state that is able to support standing waves, 
and the conjugate supercri tical flow which is unable to support standing 
waves. Subcritical flows generally have higher swirl velocities than 
the conjugate supercritical flow. In this context, the work of Squire 
gives a condition marking the interface between these two states. As in 
the work of Squire, a universal characteristic parameter was defined 
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which delineates the critical regions of the flow. This parameter, 
denoted by N, is the ratio of absolute phase velocities of long 

wavelength waves which propagate along the vortex in the axial 

direction, i.e., N=(C + + C_)/(C + -C_). Here C + and C_ are the phase 

velocities of the waves which propagate with and against the flow 
respectively. For N>1 the flow conditions are supercri tical and for 
N<1, subcri tical. 

Benjamin applied this theory to a specific vortex flow, defined by 
a constant axial velocity, U, and V=V Q r, 0<r<l and V=V Q /r, l<r<R. If 
R-*», this is just the combined vortex studied by Squire [7]. Benjamin 

found that the critical condition was the same form as Squire's. The 

precise value of the constant depends on the value of R but lies between 
1.92 when R=1 and 1.20 when R approaches infinity. Thus Benjamin, 

although starting from a different perspective, arrived at the same 
critical condition for a combined vortex as did Squire. 

A recent paper by Ito, Suematsu, and Hayase [16] examined both 

stationary and unsteady vortex breakdown in an inviscid, incompressible 
fluid. They considered the stability of a columnar vortex subjected to 
small amplitude disturbances. The disturbances considered were 
axi symmetric as well as asymmetric and either steady or unsteady. Their 
analysis produced a criterion for breakdown from the requirement for 
existence of solutions to their disturbance equations. A comparison of 
these results with those of Benjamin [8], for the case of a finite 

radius pipe containing a rigid-body rotation, gave the same criterion 
for breakdown. The important aspect of the work of Ito et al. lies in 
their interpretive criterion. Their nondimensional ization leads to the 
Rossby number as the controlling parameter. For example, in the case of 
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swirling pipe flow consisting of solid body rotation, the relevant 
scales are the axial velocity U, pipe radius r Q and constant angular 
velocity of the flow, Q. Thus the dimensionless controlling parameter 

is: R =U/r Q. 

o o 

Tsai and Widnall [17] examined a group velocity criterion which can 
be considered as a variation of the phase velocity criterion of Benjamin 
[8]. Their analyses of swirling pipe flow is more consistent with the 
view that breakdown occurs due to a wave trapping mechanism [18]. They 
assumed that the radial and axial velocity distributions could both be 
fit to exponential profiles [19]. They used a least squares fit given 
by Garg and Leibovich to calculate the dispersion relation from linear 
parallel stability theory. The group velocity associated with various 
flow profiles was then calculated. The results showed that upstream of 
breakdown the group velocity of both symmetric and asymmetric modes was 
directed downstream. Even though their criticality condition of zero 
group velocity proved to be an accurate guide for the various types of 
breakdown, they were unable to establish a relationship between vortex 
breakdown and wave trapping. 

Bossel [9] concluded that breakdown was not analogous to the 
hydraulic jump, rather it was a regular feature of the Navier-Stokes 
equations for the given flowfield. This flowfield is considered to be 
supercritical with rigid initial rotation and some axial deceleration 
near the axis. 8ossel divided the flowfield into two distinct regions: 
(1) An inner region, which could contain a stagnation point, and was 
approximated by the equations for an inviscid rotating flow; and (2), a 
viscous quasi -cylindrical region which surrounds the inner region. 
Bossel assumed the outer solution was known which produced conditions at 
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the boundary of the inviscid region that will result in breakdown. For 
a rigid rotation the inviscid equation for the stream function becomes 
linear and solutions were obtained by superposition. Results which 
resemble the configuration of a vortex breakdown were obtained. 

Hall [11] considered vortex breakdown to be analagous to the 
separation of a two dimensional boundary layer. He assumed that a 
failure of the quasi-cyl indrical approximation through large axial 
gradients signaled an impending vortex breakdown. A numerical 
experiment was performed to test the theory using experimental data 
obtained by Kirkpatrick [20]. A retardation of the flow along the axis 
was found. At this point, computations failed due to the inability of 
the iteration scheme to coverage. Hall considered this to represent the 
failure of the quasi-cyl indrical approximation. In addition, stream 
tube divergence, pressure gradient, and swirl magnitude were varied 
parametrically and effected the failure of the quasi-cyl indrical 
approximation in a manner consistent with their affect on vortex 
breakdown. Hall also found that the effect due to a change in Reynolds 
number was small. 

Linear hydrodynamic stability theory investigates only the 
amplification or decay of infinitesimally small disturbances imposed on 
the base flow. Breakdown is then assumed to be analogous to laminar- 
turbulent transition. Of course, as pointed out by Leibovich [21], 
breakdown can occur with little sign of instability and conversely a 
vortex flow may become unstable and not undergo breakdown. In the case 
of zero axial velocity the Rayleigh [22] criterion (that the square of 
the circulation should nowhere decrease as r increases) provides a 
necessary and sufficient condition for flow stability. 
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That vortex breakdown may be a result of an hydrodynamic 
instability was first proposed by Ludwig [12]. He found a stability 
boundary with respect to spiral disturbances for an inviscid rotating 
flow in an annulus. . It was proposed that this instability could be 
amplified and induce an asymmetry in the core. Ludwig's idea has not 
been widely accepted due to the difficulty of relating the geometry of 
his analyses to vortex breakdown, 

Howard and Gupta [13] have shown that the stability of the quasi 

cylindrical approximation s guaranteed if the "Richardson number" 
-1 “2 

criterion r (du r /9r) Y 0 [b(rY Q )/ar]< 1/4 is satisfied. This implies 
that the role of swirl is purely stabilizing for axisymmetric 
disturbances. In practice, nearly all approach flows turn out to be 
stable to axisymmetric disturbances [21]. 

Pedley [23] considered the stability of an almost fully developed 
viscous flow in a rotating pipe. He found that the flow became unstable 
to asymmetric disturbances for sufficiently small values of the Rossby 
number (defined in terms of the rate of rotation of the pipe, the pipe 
radius and the axial velocity of the fluid) at a critical Reynolds 

number of 82.9. 

Lessen, Singh and Paillet [14] defined a parameter, q, involving 
the ratio of the magnitude of the maximum swirl to that of the maximum 
axial velocity. This parameter completely determined the inviscid 

stability characteristics of the flow defined by the equations 

-r 2 -r 2 

V = ( q/r ) ( 1-e ), U - e . Thus a wake or jet-like axial velocity 

profile does not affect stability. They found that the flow was stable 

to all disturbances for q>1.5 and unstable to nonaxi symmetric 
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disturbances for smaller values of q. Thus, from the point of view of 
stability, the role of swirl is stabilizing with respect to axisymmetric 
disturbances, and destabilizing with respect to asymmetric disturbances 
over a range of q. 

After the initial observations of vortex breakdown by Peckham and 
Atkinson [3], experimentalists began studying vortex breakdown in a more 
controllable setting. Harvey [24] performed experiments in a long tube, 
imparting a swirl velocity on the fluid as it entered. This was done 
using a set of adjustable vanes mounted in the inlet section. Harvey 
concluded that for low swirl velocities the classical vortex was 
obtained but as the swirl was increased a breakdown bubble formed. He 
also concluded that the breakdown was due to a critical state phenomena 
and not a hydrodynamic instability since the flow reverted to a normal 
swirling flow downstream of the breakdown bubble. Instabilities usually 
result in increasingly large amplitude velocity fluctuations ending in a 
turbulent flow. 

Sarpkaya [25] described experiments in swirling flows in a 
diverging cylindrical tube. He observed three types of breakdown; 
double helix, spiral, and axisymmetric. The type of breakdown that 
occurred depended on a combination of Reynolds number (based on tube 
diameter and mean axial velocity) and circulation. For 1000<Re<2000 the 
spiral or double helix breakdown was observed. Axisymmetric breakdown 
was found to develop from the double helix or spiral form, or as an 
axisymmetric swelling of the core. For high Reynolds numbers and 
circulation the axisymmetric type breakdown occurred as a swelling of 
the vortex core. He also noted that the axisymmetric breakdown 
responded to changes in upstream and downstream flow conditions in a 
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manner analagous to hydraulic jump in an open channel flow. In a later 
paper, Sarpkaya [26] concluded that Benjamin's finite transition concept 
was in complete agreement with experimental results in the region where 
axisymmetric breakdown occurred (high swirl and Re). The spiral 
breakdown appeared to be a consequence of the instability of the flow 
due to asymmetric disturbances. He concluded that the overall mechanism 
for vortex breakdown might encompass finite transition and hydrodynamic 
instability theories, each applicable in a specific region. 

Faler and Leibovich [27] have mapped the internal structure of an 
axisymmetric type vortex breakdown using a laser-Doppler anemometer. 
They found that the interior of the bubble, which contained a two celled 
structure, was dominated by low frequency periodic velocity fluctua- 
tions. The magnitude of these fluctuations was greatest in the rear 
portion of the bubble. In addition, four stagnation points existed on 
the axis. 

The affect of an adverse pressure gradient on vortex breakdown has 
been examined by Sarpkaya [28] and more recently by Delray et al. 
[29]. Delray found experimental limits for vortex breakdown as a 
function of adverse pressure gradient and vortex strength. Pressure 
measurements showed considerable pressure increase within the core for 
small variations outside the core. In general, as the adverse pressure 
gradient increased, the swirl required to induce vortex breakdown 
decreased. 

Numerical solutions for vortex bursting have been reported by at 
least seven previous investigators [30,31,32,33,34,35,36]. In all 
cases, the flows were assumed to be incompressible and were restricted 


10 


to studies of laminar, axisymmetric systems. The solutions by Kraus 
et al. and Hafez et al. were time dependent while the others were steady 
state solutions. 

Kopecky and Torrance [31] considered axisymmetric swirling flow 
through a cylindrical tube. The fluid entered the tube with an 
exponential swirl velocity. This distribution behaved as a solid body 
near the axis and as a potential vortex away from the axis, representing 
a solution to the Navier-Stokes equations for the limiting case of 
Reynolds number approaching infinity. A parametric study was performed 
with Reynolds numbers (based on tube radius and constant axial velocity) 
ranging from 50 to 500 and swirl ratios from 0.4 to 10. The development 
of a recirculation zone was demonstrated as the swirl was increased for 
fixed Reynolds number and core diameter. Similar results were obtained 
when the core diameter and swirl were fixed while the Reynolds number 
was increased. In all cases the breakdown appeared to form at the inlet 
station. With a grid spacing of 0.25 in the streamwise direction the 
major portion of the breakdown was contained within about four grid 
points. While this grid seems excessively coarse, Kopecky and Torrance 
reported that doubling the number of grid points in the streamwise 
direction produced similar results. 

Grabowski and Berger [33] solved the steady axisymmetric Navier- 
Stokes equations for a free vortex approximated by a two parameter 
family of assumed inflow distributions. The inflows were the polynomial 
profiles given by Mager [37] in his integral analysis, imbedded in an 
irrotational flow. The equations of motion were written in terms of 
stretched coordinates in the radial and axial directions. At inflow, a 
parameter, alpha, could be varied to allow for jet like or wake like 
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axial profiles. An "artificial compressibility" technique was used to 
solve the equations of motion. Solutions were obtained which exhibited 
many of the characteristics of vortex breakdown for Reynolds numbers up 
to 200. These solutions were obtained with upstream conditions that 
were, in many cases, subcritical. Their results appear to refute the 
finite transition theory of Benjamin [8] which required the flow 
upstream of breakdown to be supercri tical . 

Narin [32] investigated the occurrence of vortex breakdown for 
three different flow configurations: (1) a straight tube, (2) a step 
tube, and (3) confined jet mixing. This appears to be the only work 
investigating the breakdown of a confined jet, which consists of a 
swirling jet discharging into a coaxial nonrotating surrounding 
stream. For this configuration, the resulting flow field depended on 

the radius of the enclosing tube, jet velocity and swirl ratio, and on 
the velocity of the surrounding stream. In general, increasing Reynolds 
number and swirl ratio enhanced the severity of the vortex breakdown. 

Benay [35] has also simulated vortex breakdown by a numerical 
solution of the laminar axi symmetric Navier-Stokes equations. At 
inflow, an exponential circumferential velocity distribution was imposed 
in a parametric study to determine the effect of vortex core radius, 
Reynolds number (based on tube radius and free stream axial velocity), 
and tangential and axial velocity on vortex breakdown. In general, an 
increase in the Reynolds number or swirl ratio resulted in a more 
pronounced recirculation zone. 

Excellent review articles summarizing vortex breakdown research 
have been published by Hall [38] and Leibovich [21,39]. Since relevant 
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numerical work only began in the mid 1970's, only Leibovich's survey 
includes those investigations. He noted that the numerical solutions do 
not show a two celled structure within the breakdown bubble as revealed 
in physical experiments by Faler and leibovich [27]. In addition, Faler 
and Leibovich measured four stagnation points along the axis, whereas 
numerical experiments have shown only two. However, this structure was 
later claimed to have been computed by Kraus, Shi, and Hartwich [34] by 
studying the flow in a time dependent manner. An examination of their 
computed streamlines reveals that the bubble has lifted off the axis for 
some of the time levels shown. At these time levels no stagnation 

points are present along the axis. Leibovich is also critical of the 
fact that the numerical solutions “contain strong axial gradients right 
up to the initial axial station." He suggests that the bubble may pass 
through the initial stations if the inlet boundary conditions were 
relaxed. In addition, axisymmetric numerical solutions show bubbles 
that increase in size as swirl is increased, a result that is not 
consistent with experimental observations. Furthermore, axisymmetric 
numerical simulations are unable to predict spiral type breakdowns, 
which are a common occurrence in experiments [6,25]. Thus, Leibovich 
has concluded that "the assumption of steady axisymmetric motion may not 
be adequate to compute all the detailed structure of vortex breakdown." 

The purpose of this work was to study numerically the spatial and 
temporal evolution of a class of vortical structures. Wing tip vortices 
are of specific interest, but the influence of wing geometry on the tip 
vortices was beyond the scope of this analysis. A numerical solution 
has been chosen because closed form analytic solutions of the equations 
.of motion are unlikely to be found without assuming overly restrictive. 
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simplified flow fields. An experimental approach was also deemed 
impractical because of the difficulty in making accurate measurements in 
regions where the velocity gradients are extreme, which is often the 
situation during vortex breakdown. Thus, the unsteady, three- 
dimensional, laminar Navier-Stokes equations have been integrated 
numerically to study the parameters affecting the evolution and possible 
breakdown of an isolated wing tip vortex. 

Authors of the numerical studies cited previously report breakdown 
at or immediately downstream of the inflow boundary. It will be shown 
that these previous computational results can be re-evaluated in terms 
of a single parameter which identifies the cause of the breakdown at 
inflow. Through an examination of several previous studies 

concentrating on standing wave analyses, it is shown that the 

controlling parameter is the Rossby number. A means of avoiding the 

problem of breakdown at inflow is suggested, and the numerical analyses 
is then performed with the Rossby number and Reynolds number as the 
nondimensional parameters. It is important to emphasize the fact that 
the algorithm was not restricted by an axisymmetry requirement. This 
allowed for the existence of asymmetric disturbances which may be 
important in the breakdown process. Previous numerical simulations have 
precluded this possibility. The numerical algorithm which has been used 
is an implementation of the "veloci ty-vortici ty" formulation of Gatski, 
Grosch and Rose [39]. The pressure, although not a variable in the 
formulation, was computed. 

Contour plots of pressure, Bernoulli's constant, axial velocity and 
vorticity are displayed as a part of this investigation. Particle 
traces, vortex lines and velocity vector plots have been obtained using 
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the PL0T3D, a three-dimensional color graphics program implemented on an 
Iris workstation at NASA Langley Research Center. In addition, the 
rates of change of energy and enstrophy were computed and plotted as a 
function of axial location. Conclusions have been drawn which may be 
useful in the interpretation, as well as the modification and control, 
of wing tip vortices. Finally, the internal structure of the vortex 
breakdown bubble is discussed. 


15 



THE ROSSBY NUMBER - BREAKDOWN CRITERION 


The Rossby number is an important control parameter in the study of 
large scale atmospheric and oceanic motions. It is a measure of the 
relative importance of the Coriolis and inertial forces on the fluid 

motion. The Coriolis force is due to the rotation of fluid and is 

directed perpendicular to the axis of rotation. In the study of 

geophysical fluid dynamics, fluid rotation is generally considered to be 
of the rigid body type. However, the Coriolis forces can be important 
for a variety of circumferential velocity distributions associated with 
other flows, such as those occurring in wing tip and leading edge 
vortices. When significant, the Coriolis acceleration represents a 
restoring force, providing a mechanism for the creation of waves (in the 
absence of sufficient damping). It tends to restore fluid particles 
displaced laterally from their equilibrium positions. However, the 
restoring force can cause the fluid particles to overshoot their 
original locations, setting up an oscillatory motion. Under some 

conditions the fluid can sustain these oscillations, and in the case of 
vortex flows, these wavelike fluctuations can then propagate along the 
axis of the vortex. Waves of this type are known as inertial waves 
[40]... The intent of this chapter is to show how the vortex breakdown 
phenomenon can be characterized in terms of the ability of a base flow 
to support these waves. This effect can be described, in terms of the 


Rossby number, and can be justified by using the theoretical, experi- 
mental and computational literature discussed in the introduction. 

Throughout the remainder of this chapter a cylindrical polar 
coordinate system, (r, 9, x), with corresponding velocity components, w 
in the radial (r) direction, v in the circumferential (9) direction, and 
u in the axial (x) direction are employed. 

The Rossby number can be developed naturally from the vorticity 
transport equation and is defined as the ratio of the inertial forces to 
the Coriolis forces as 

Ro * -4- (2.1) 

r Q 

★ 

where u is a representative velocity magnitude, r a characteri stic 

length, and Q, a characteristic rate of rotation of the flow. For the 

★ 

flows considered in this study, r is taken to be the vortex core 

diameter, defined as the radius of maximum swirl velocity. The 

reference velocity u is taken as the axial velocity at the core radius 
★ 

(r ). Wing tip vortices are characterized by a solid body type rotation 
hear the axis, and this is taken to be the' characteristic rate of 
rotation, Q. 

As discussed earlier. Squire [7], Benjamin [8] and Ito et al. [16] 
were able to derive a criterion for the existence of standing waves on a 
vortical base flow. Squire [7] and 8enjamin [8] formulated this 
criterion in terms of characteristic circumferential and axial 

velocities. Here it will be shown that their criteria can be 
reinterpreted in terms of a Rossby number. 
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The most important case considered by Squire, in terms of a model 
for vortex breakdown, was an exponential form for the circumferential 
velocity profile, given by: 



with V Q used as a scaling parameter. The axial velocity was considered 
to be a constant, i. e. u = U. Recall from the introduction that the 
existence of neutrally stable standing waves occurred when 

V/U = 1.0 (2.3) 

m 

where V m is the maximum swirl velocity. In addition, it should be noted 
that the maximum swirl velocity using Squire's velocity profile is 
V m = 0.638 V 0 , at r* = 1.12. 

Consistent with our previous definitions, the reference length, r* , 
is given as: 

r*» 1.12. (2.4) 

The reference velocity is the constant axial velocity, U. The charac- 
teristic rate of rotation, Q, is given as: 

Q = lim (v/r) * V . (2.5) 

r*o 0 

Hence the Rossby number is then given as: 

Ro = 0.57. (2.6) 


18 


The combined vortex considered by 3enjamin [3] is given by 


v = V r 
o 

V - V 0 /r 
u = U 

where V Q is a constant. 

The critical condition, for R=1 and R 

V Q /U = 1 .92 
V Q /U = 1.20 

The case R*®, corresponds to the 
[7]. The case R=1 corresponds to 
and was also studied by I to et al. 
taken to be the distance at which 
tional flow are matched. Thus, r* 
V Q , These results can be expressed 


0<r< 1 

l<r<R (2.7) 

all r 

* « is given as 
for R = 1 

( 2 . 8 ) 

for R .-*> » 

combined vortex studied by Squire 
a solid body rotation within a tube 
L 16 J . The characteristic radius is 
the solid body rotation and irrota- 
is equal to unity and 2 is equal to 
in terms of a Rossby number as 


Ro * 0.52 for R * 1 

Ro * 0.83 for R * ® 


(2.9) 


From the above analyses, it appears that a criterion based on the Rossby 
number can be used to delineate the critical regions of the vortex flow 
but the critical value depends upon the type of vortex flow. This 
criterion was then used as a basis for examining a variety of vortex 
flows. Previous computational and experimental work has been examined 
for both confined and unconfined flows to determine the range of 
applicability of this Rossby number parameter. 
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The computational studies of Kopecky and Torrance [31], Grabowski 
and Berger [33], Benay [35] and Hafez et at. [36] have been 
reinterpreted in terms of the Rossby number. The circumferential 
velocity profiles used by Kopecky and Torrance and Benay are of 
exponential form, similar to Eq. (2.2). Grabowski and Berger [33] and 
Hafez et al. [36] express the circumferential velocity in terms of a 
polynomial. Both profile types asymptote to solid body rotation near 
the axis. The exponential profile asymptotes to an irrotational flow in 
the far field. The polynomial profile is exactly irrotational outside a 
specific core radius. 

The circumferential velocity .profile of Grabowski and Berger [33] 
and Hafez et al. [36], in nondimensional form, is expressed as 


v = Vr( 2 - r 2 ) 0<r<l 

v 3 V/r 1 <r<R 


The axial velocity profile is given as 

2 

u 3 a + (l-a)r(6-8r+3r ) 0<r<l 

u * 1 l<r<R 


( 2 . 10 ) 


( 2 . 11 ) 


where a is an adjustable parameter to allow for jet-like or wake-like 
profiles. The circumferential velocity is a maximum at r = /2/3, and is 
equal to 1.088 V. The characteristic rate of rotation, Q, is given as 

Q = lim (v/r) 3 2 V U /6 (2.12) 

09 

r-*-o 
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where 5 is the dimensional core radius at the inflow plane and 1) the 

co 

dimensional free stream axial velocity. The radius of maximum swirl 
velocity is given as 

r* = /m 5. (2.13) 

The characteristic axial velocity becomes 

u = {a + (1-a) 8/273 (1- /2/3)} U (2.14) 

00 

For the case of a 5 l,i.e. u = , the Rossby number becomes ■ 

Ro = .6 12/ V (2.15) 

where V is a parameter describing the circumferential velocity 
profile. The Reynolds number, based on the free stream axial velocity, 
U^, and the character!' Stic radius, r* , becomes 

Re a /273 U 5/ v (2.16) 

09 

where 5/v is the form of the Reynolds number employed by Grabowski 
and Berger [33] and Hafez et al. [36], 

In a similar manner, expressions for the Rossby number and Reynolds 
number can be extracted from the works of Kopecky & Torrance [31] and 
Benay [35]. Using the notation employed by these authors, the following 
results are obtained for u = constant 

Ro (2.17) 

1.12 /T r . 

o 

Re = — (2.18) 

/T (U r /v). 

® o 
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The results of the previous works are summarized in Fig. 1.1 in 
terms of the Rossby and Reynolds numbers. Open symbols denote breakdown 
and closed symbols denote no breakdown . That figure and Fig. 1.2 have 
employed an axial velocity criterion to define breakdown. For confined 
flows and numerical solutions, vortex breakdown is defined as a flow 
which produces stagnation of the axial velocity component. For 
unconfined flows, breakdown is considered to be a rapid expansion of the 


vortex core 

accompanied 

by a 

sudden 

deceleration of 

the 

axial 

ve 1 oc i ty . A 

1 imi t 

1 ine 

exists 

that 

separates regions 

of 

vortex 

breakdown from 

regions 

that 

experience no 

breakdown. 




The computational results show a Reynolds number dependence in the 
low Reynolds number range. From Fig. 2.1, it can be seen that for 
Reynolds numbers above 100, viscous effects appear to be negligible and 
inviscid theory can be expected to give good results. The dashed line 
represents the inviscid standing wave theory developed by Squire [7] for 
an exponentially varying circumferential velocity profile. 

The experimental results appearing in Fig. 2.1 are all at higher 
Reynolds numbers than the numerical simulations. Garg and Leibovich 
[19] and Uchida et al. [41] made LDV measurements just upstream of 
breakdown in a tube and vane apparatus. In the case of Garg and 
Leibovich, a least squares fit of the data was used to obtain 
exponential profiles for both the circumferential and axial velocity 
distributions. These results were easily translated to the Rossby 
number and Reynolds number previously defined. The bubble form of 
breakdown occurred at a lower Rossby number (-0.57) than the spiral 
form (~ 0.63). From the available data, it appears that the Rossby 
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Fig. 2.1 Correlation of Rossby number for vortex breakdown (shaded 
symbols) with Reynolds number. Wing tip vortex case. 


O Owen & Peake 
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Fig. 2.2 Correlation of Rossby number for vortex breakdown (shaded 
symbols) with Reynolds number. Leading edge vortex case. 


number at inflow, in the case of the bubble form of breakdown, was 
already below the level of the Rossby number at which the spiral form of 
breakdown occurred. The Rossby number for the Uchida et al. work, 
obtained from plots of the axial and circumferential velocity profiles, 
equaled 0.64 for the bubble form of breakdown. Note the excellent 

agreement between these confined experimental flows and the inviscid 
standing wave analyses of Squire. 

Singh and Uberoi [42] measured the velocity distribution of a wing 
tip vortex at several axial stations along the vortex core. A laminar 
flow wing was used to generate the vortex. Their measurements 

provided enough information to obtain an estimate of the Rossby number, 
(~ 0,60) and Reynolds number, (~ 13000) just upstream of a region in 
which the axial velocity decreases rapidly to 0.3 U a , suggesting vortex 
breakdown . 

8ased on Fig. 2.1 the critical Rossby number for the symmetric form 
of breakdown for the trailing wing tip class of vortices is 
approximately 0.60. The computational results indicate that for 

Reynolds numbers below 100, the value of the critical Rossby number is 
decreased, undoubtedly due to the increased damping effects of viscosity 
on the wave motions. 

Figure 2.2 displays the Rossby number-Reynolds number relationship 

for leading edge vortices. The experimental data were obtained from 

reports by Owen and Peake [43], Anders [44] and Verhaagen and Krui shrink 
[45]. Once again, open symbols denote no breakdown and closed symbols 
denote breakdown. 
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Anders [44] made LDV measurements of a leading edge vortex over a 
delta wing at two different angles of attack; 19.3* and 28.9*. The 
vortex produced by the wing at the lesser angle of attack did not break 
down. The vortex produced at the angle of attack of 28.9“ degrees broke 
down above the delta wing. The Rossby number is computed at the same 
distance from the apex of the wing for both cases. 

Blowing can stabilize a leading edge vortex and thus change the 
breakdown criterion. Owen and Peake [43] introduced core blowing into 
vortices shed from delta wings at high angles of attack to study its 
effect on vortex breakdown. The symbols in Fig. 2.2 representing these 
data are variations based on a blowing coefficient parameter, C , at 
fixed streamwise locations z/c=3 and z/c=4, where c is the chord length 
of the delta wing. Owen and Peake state that breakdown occurs 
for the case C = 0.0, while for the cases C = 0.05 and C * 0.12 the 
flow is stabilized and no breakdown occurs. 

Although the data are sparse and the evaluation of the Rossby 
number approximate, one may conclude that vortex breakdown for leading 
edge vortices occurs at a higher Rossby number than for trailing wing 
tip vortices. This may be due to the fact that the swirl velocity 
profiles are of a different type. Far downstream, the flow outside the 
core of a trailing wing tip vortex is nearly irrotational . For a 
leading edge vortex, the flow at the edge of the core is rotational and 
nearly inviscid. In addition, the leading edge vortex contains a narrow 
viscous subcore where the radial gradients of the circumferential 
velocity are extremely large. In contrast, a significant region of 
fluid near the center of a wing tip vortex tends to rotate like a rigid 
body. 
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Mo standing wave analyses of velocity profiles applicable to 
leading edge vortices were found during this investigation. That 
analysis could predict an analytic Rossby number breakdown criterion for 
leading edge vortices and would bridge the two cases.’ However, based on 
experimental results, the critical Rossby number should be near unity. 

In summary, the theoretical analyses of Squire [7], Benjamin [8] 
and Ito et al. [16] have been reinterpreted to enable the identification 
of a criterion that predicts the existence of axisymmetric standing 
waves based on a Rossby number. An exponential representa tion of the 
circumferential velocity profile, which most closely models experimental 
flows, yields a critical Rossby number of 0.57. This value is shown as 
a dashed line on Fig. 2. The experimental data of Garg and Leibovich 
[19], interpreted in terms of a Rossby number, shows that the spiral 
form of breakdown occurs when the local Rossby number falls to 
approximately 0.63. From the available data, the local Rossby number 
was initially below 0.63 for the cases involving the bubble form of 
breakdown. Numerical experiments reveal that the critical value of the 
Rossby number for the bubble form of breakdown' becomes independent of 
Reynolds number above Reynolds numbers of 100. At lower values of the 
Reynolds number, a lower Rossby number is required to initiate 
breakdown. 
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NUMERICAL FORMULATION 


The governing equations describing incompressible , isothermal flow 
of a Newtonian fluid are given by: 



7*u 

3 0 


(3.1) 


Du 

-1 

2 - 





7P + v7 u 

(3.2) 


Dt~ 

P 



subject to 

u 

= 9 

on boundary 8. 

(3.3) 

Here, Eqs. 

(3.1) and 

(3.2) 

represent the continuity 

equation and 


the Navier-Stokes equation's respectively, each valid over a domain D. 
Equation (3.3) is a statement of specified velocity boundary condition 
to be satisfied on the boundary, 8, of the domain. Higher order 
equations involving the vorticity, id, are given by: 


7«u = 0 


7xu = <D 

a U)*7u + v7^ (D 


7 3 0 

with a corresponding boundary condition: 

(D 3 7 X U. 


(3.4) 

(3.5) 

(3.6) 

(3.7) 


7 Q 


on 8. 


Equation (3.4) represents the definition of vorticity, Eq. (3.5) the 
vorticity transport equation and Eq. (3.6) the solenoidal condition on 
the vorticity vector. Equation (3.6) is an identity obtained by taking 
the divergence of each side of Eq. (3.4). The numerical scheme used to 
solve these equations represents an implementation of a method described 
by Gatski, Grosch and Rose [39]. The scheme is second order accurate in 
time and space for the vorticity and velocity variables, and is 
independent of the local cell Reynolds number. A major advantage of 
this formulation is that boundary conditions for the pressure are not 
needed to advance the solution in time. (The difficulty in specifying 
the pressure boundary conditions accurately has been discussed by 
Orszag, et al . [46].) 

Several major aspects of this algorithm can be identified. 
Equations (3.1) and (3.4) form the basis for the solution to the 
velocity vector field when given the vorticity vector field at any time 
level n, along with the velocity boundary conditions. Equation (3.5) is 
utilized to advance the vorticity field from time level n to time level 
n+1. Here, the boundary condition, Eq. (3.7),- is needed to produce a 
unique solution to the higher order system of equations which are 
solutions to Eq. (3.1) and Eq. (3.2) subject to the boundary condition 
prescribed by Eq. (3.3). The solution of Eq. (3.6) is used periodically 
during the time evolution of the flow. field to ensure that the vorticity 
vector remains divergence free. That is, the divergence free 

requirement is tested at each time level and if it fails to meet its 
tolerance, Eq. (3.6) is employed. 

For reference purposes, a brief description of the computational 
sequence follows. The physical domain is first divided into a 
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computational region of rectangular cells IMAX x JMAX x KMAX. A typical 
cell is sketched in Figure 3.1. Velocities, defined at the centers of 
faces, are the average of box variables defined at the vertices of 
cells. The numbering scheme (1-8) for the box variables is shown in 

Figure 3.1. Vorticity is also defined at the center of each cell 
face. The velocity and vorticity variables thus represent average 
values over a cell face. 

Beginning with an assumed vorticity distribution, the velocities, 
utilizing Eqs. (3.1) and (3.4) are computed at time level n. The 
v-orticity is then advanced to level n+1 using the velocities at time 
level n and vorticity boundary conditions determined by the velocity 
components on each boundary. This vorticity is subsequently projected 
into a new vector space satisfying the requirement that the divergence 
of vorticity be zero. Next, velocities are updated to time level n+1 
usirrg the divergence free vorticity at n . and appropriate velocity 
boundary conditions. The vorticity is then recomputed at level n+1 
using the updated velocity field. Repetition of the above process 
yields a second order accurate solution at any subsequent time [39]. 

Numerically, it is required that one component of the velocity 
vector be specified on each boundary cell face. All three components of 
the vorticity vector are specified on each boundary cell face. 

Velocity Equations 

Assuming the vorticity at time level n is known, the velocity 
components at level n can be computed through a numerical solution of 
Eqs. (3.1) and (3.4) where u satisfies the compatibility condition, Eq. 
(3.6). Fix and Rose [47] have shown that this compatibility condition 
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Fig. 3.1 Representative computational cell showing the location of 
the velocity and vorticity variables. 
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is a necessary condition for the overdetermined system, Eq. (3.1) and 
Eq. (3.4) to have a solution. 

A finite volume approximation to the continuity equation is 
obtained by applying the divergence theorem to the integral form of the 
continuity equation. That approach results in the requirement: 

/v*udV * /u*nda (3.3) 


The integral on the right hand side represents the net flux of mass 
into any arbitrary fixed volume. It is convenient to define an operator 


V u • as 


V 

h 


•u 



/ u *n da. 


(3.9) 


Using the trapezoidal rule, the continuity requirement can be expressed 
as 

7 h *u = Z 6.u. i = 1,2,3 (3.10) 


(where 6. is the standard difference operator, 5 u = (u , -u , )/&x ) 
1 i i i+V2. i - Hi i 

so that 7 h «u » 7*u + 0 (h 2 ) ' (3.11) 


The discretized form of the definition of vorticity results from an 
application of Stokes theorem. Consider a two-sided surface in three 
dimensions having a closed surface C as its boundary. The circulation 
of the velocity u around C is equal to the flux of the curl of u over S, 
i .e. , 


/ (7xu)*ndS - j u*dC 
s J c 


(3.12) 
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where n is a unit normal to the surface S. We can define a second 


opera tor 7^ x u as 

7 h x u 'tit ^ u-aC ua3) 

Using the trapezoidal rule, the discrete form of the curl operator is 
given as 

(7 x u) 3 5. u. - 6.u. (3.14) 

h k i J J i 

so that V^xu=Vxu+ O(h^) (3.15) 

The preceding approximations are valid over the entire 
computational domain. Additional equations are needed for boundary 
cells. For the numerical solution, a Dirichlet condition must be 
applied to any one component of velocity on all boundaries. This 
condition is expressed as 

Jg u*n dS = 0 (3.16) 

The numerical problem requires the solution of Eq. (3.10) on each 
cell, Eq. (3.14) on each cell face, and Eq. (3.16) on all boundary cell 
faces. This can be solved by an iterative scheme due to Kaczmarz 
[48]. The resulting solution represents a least squares approximation 
to the system Ax = F, where A is an nxm matrix and x and F are m and n 
dimensional vectors respectively. When applying this scheme to the 
system of equations given by Eq. (3.10) and (3.14), each equation is 
relaxed independently. Therefore, n = 1, and m is equal to the number 
of unknowns in a cell. The scheme is derived in the following 
paragraph. 
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If R is a residual and (k) ,is a counter, then we can write: 

A^u^- F. = R..^ and u^ k+ ^ = u^ + v^j, where v is a correction term. 

- F. where i represents the specific equation to be 


R. 

i 


1 1. i u _ _ r t, j. i > 

\ f.-r 1/ _ ~ 


= A. «u' 


relaxed. 


. A, (i (kl ♦ : (k) ) - F, 


1 


1 


choose 


= R. (k) + A i -v (k) 


Rj (k+1) . 0 


0 = R . ^ k ^ + A. v (k) ; define v (k) * a! w (k) . 


Then 


A i A. w (k) = -R. (k) 


w (k) * -<A A. 1 )' 1 R< k) 


v (k) = -a! (A. A .)" 1 Rj k) 


-(k+i) 3 -(k) (A.A i f 1 R| k) 


(3.17) 


An acceleration parameter y has been introduced in Eq.(3.17). 
Expressions for (AA^) appear below for each of the required equations. 

■I 


continui ty 


x vorticity 


(A A T ) 


(A A T ) 


-1 


1 

7~"7“ 

3(l+a +r) 


1 


4(1+kt 4 ) 


(3.18) 


(3.19) 
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y-vortici ty 

_ _t 1 

(A A 1 ) = 

4(l+fT) 

(3.20) 

z-vortici ty 

(A A T . —L_ 

4 ( 1+a ) 

(3.21) 


where < = Ay/Az 

(3.22) 


6 = Ax /Ay 

(3.23) 


and a = Ax/A z 

(3.24) 

Dirichlet type boundary condition 



- -T 

(A A 1 ) = 1/4. 

(3.25) 

To implement 

the Neumann condition, du/dx = 0 on a cell 

face, the 


continuity equation can be used to write 

3v/3y + 3w/3z = 0. Then 

. . T 'I 2 

(A A 1 ) * .25/(1+* ) (3.26) 

In summary, a projection method due to Kaczmarz has been 
implemented to solve the overdetermined system given by Eqs. (3.1), 
(3.3) and (3.4). Tanabe [49] ha-s shown that the method will converge 
for any system of linear equations with nonzero rows even if the system 
is singular. 


Vorticity Equations 

The discretized form of the vorticity transport equation is 
obtained by expressing the' vortici ty within a cell in terms of a set of 
basis functions. These functions are integrated over time and space 
resulting in expressions valid on cell faces. Appropriate combinations 
of these expressions result in the discretized form of the equations for 
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the transport of the vorticity. This procedure is detailed by Mclnville 
[50] . 


The basis set employed by Gatski, Grosch, and Rose [51] for the 
two-dimensional formulation is a solution to the one 'dimensional form of 
the vorticity transport equation. These solutions are of the form 


cjj(x,t;a) = exp [ax-0(a)t] 
where 0(a) = a(u - av) 


(3.27) 

(3.28) 


The three-dimensional form of the vorticity transport equation 
contains the vortex stretching term u»7u which requires special 
consideration. The basis set no longer represents solutions to the 
three dimensional transport equation. Thus, Gatski, Grosch and Rose 
[39] employed the transformation 


CO 




(3.29) 


where co is the vorticity and C is a transformed vorticity. Matrix [B n ] 
is a 3x3 array which relates each component of Z to the three components 
of C. When applied to the vorticity transport equation, this 
transformation eliminates the vortex stretching term. The resulting 
transformed equation is 

§|- + u«7C = vT 2 C (3.30) 

Details of the transforma tion appear in Appendix A. This form of the 
transport equation has solutions of the type (3.27) and is the equation 
to be discretized. Note that at time level n the transformed vorticity 
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equals the physical vorticity. Hence, "true" vorticity and 
"transformed" vorticity are equated at the beginning of each time step. 

The discretized form of the transport equation has been derived in 
Gatski, Grosch and Rose [39]. It takes the form. 


where. 


{6 + [(n u)6 +(tiv)6 +(uw)5]}c n 
t xx y y z z 


= v (5 $ + 5 4 + 5 £) 

x y z 


(n - hq5)0=6C 

X XXX x 


(u -hq6)<l>=5C 

y y H y y y w 


( “z - h zVz> 1 * 6 z' 


u c = H C - h p 6 0 
t X XXX 


*t c 


n C * h p 6 4; 

y y H y y 


“t' 


“z 5 - h z p z 6 z 5 


K. . 5C, . SC. . 

sr 1 +ir- s +sr k 


7 K i : 3C 2 : 3C 3- 

* + 3 T J *~ k 


ac ac ac 
- i : 2 : 3 ; 

5 * - — i + - — j + k 

az az az 


(3.31) 


(3.32 a,b,c) 


(3.33 a,b,c) 


(3.34) 

(3.35) 

(3.36) 
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and 


(3.37) 


q(9 z ) = coth ( Q j) - 9' 1 ; l = 1,2,3 


p ( V = p ( V /9 a 


(3.37a) 


9„ = 




jf ; * * 1,2,3 


(3.38) 


The finite difference operators 5 and n are defined as follows. 


s n - s n 

5 S n = l /2 .J.k 

x i,j,k Ax 


(3.39) 


5 s". 

t i.j.k 


s n+ ^2 _ s n ~ \ 
l.j.k i , j , k 

At 


(3.40) 


s n + ^n 

s n = 21^-lllLl i “ 1/ 2»3> k 

x i.j.k 2 


(3.41) 


c n 

i.j.k 


s" + ^ * s"‘ 1/2 

ijk ijk 


(3.42) 


The algorithm governing- the time advance of vorticity from time 
level n to level n+1 utilizes the vorticity at time level n+ V 2 which is 
given by 


- n+ 

0) 


B n (t “ O 

n+ V 2 -z n+ V2 


(3.43) 


By definition 



(u + t 6 )C n 


(3.44) 
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Thus 

Now consider 
n+1. 


the 


3 (t - O 
- n+ V 2 = e n+ V 2 

vorticity at time level 


(n t + x6 t )C n (3.45) 

n+ V 2 centered around level 


. -8 (t - t 1 1 

- n+ V? n+1 n+1 n+ V? )- n+ V? 
co 4 * e c C 


(3.46) 


By definition 



5 \ )c 


n+1 


(3.47) 


Thus , 


- n+ V 2 
u 


= e 


-B (t - t . ) 
n+1 n+1 n+ V 2 


(n - 

t 


- n+1 
6 )C 
t 


(3.48) 


By equating Eqs. (3.45) and (3.48) the condition that the vorticity at 
time level n+ Vj? be continuous is imposed. This resul ts in a relation- 
ship between the vorticity at time level n+1 and level n. This 
condition is given as: 

, 8 t" 8 t 

, , “ n+1 n+1 r n , , - n_ 

(u - x 6 ) c 3 e [e (u + t6 ) C ] (3.49) 

U w U L 


where 



(3.50) 


*=t„ +V2 -t n (3.51) 

Equation (3.49) governs the advance of vorticity from time level n 
to level n+1. If it is assumed that the vorticity field is known at 
time level n, then the right hand side of Eq. (3.49) can be calculated 
explicitly. The time average on the right hand side is expanded by 
using Eqs. (3.33 a,b,c). The time difference is expanded by using the 
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transport equation ( Eq. (3.31)). These expansions remove the n+ V 2 
time level from the right hand side of Eq. (3.49). The expanded forms 
of Eq. (3.49) using Eq. (3.31) and Eq. (3.33a) (3.33b) or (3.33c) 

respectively, become: 

, 8 x 

, _ - n+1 n+1 n r - n . -n -n , -n_ 

(p - -c 5 )C = e e Cp C * x [ji u6 C + p v6 C + p w6 C ] 
tt x xxyyzz 

— _ 2 • n 

+ vt(5 0 + 5 4 + 5 £) *h p 6 0 ] (3.52) 

x y z xxx 


( V T 'V‘ 


n+1 


8 t" B x 

n +1 n r -n . , -n -n -n, 

e e [11 C * x[ p u 6 C + p v 6 C + p w 6 C J 
y xx y y z z 


+ vx(5 x 


(1) + 6 dj + 
v 



P .A*"] 


(3.53) 


(p t ~x'S t )C 0+1 * e n+1 e n [p 2 C°- x [p x u5 x C n + P y v6 y C n + ^w6 z C n ] 


+ vx (5 0 + 5 4+5 5 ) - p'5^ n ] (3.54) 

x y z z z Z 

■The only unknowns on the right hand side in these equations are the 

diffusion terms 6 ^ 4 , 6^4, and 5^1. Explicit expressions for these 

terms are obtained from a set of equations, derived from the identity 

- n- V? - n - n — n~ V? 

C c = p C • x 6 C where C is known from the previous time 

t t 

step. This identity is expanded using Eq. (3.31) and Eqs. (3.33 a,b,c). 
The three resulting sets of equations become 


^ =p c n +x(p u 6 C n +P v 6 c n +P c n )“(h ^p +vx )6 4~vx5 4-vx5 c, (3.55) 
x xxyyzz xx x y z 
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n u ^ n +T (p, U 5 c n+ u v6 C n+ u w6 c n )~v-c6 <&- { p +vt)<5 <Jr-vx6 £ (3.56) 

y x x y y z z x yy y z 


r n ^ a c n +x(p. u6 C f1+ M’ v6 C +u w6 C )-vc6 $-vt 6 4r-(h“p +vx)5 Z, (3.57) 
z x x y y z .z x y z z z 

Note that the i component of the above equations represents one "set," 

consisting of 3 scalar equations with unknowns 5^, 5^^ and 6^. 

Similarly for the j and k components. 

The details of the solution of these equations for 5 0 , 5 ii>, 6 5 

x y z 

by utilizing Eqs. (3.32 a,b,c), appears in Appendix 8. The advance to 

time level n+1 is implicit. Equations (3.52), (3.53) and (3.54) are 

reduced to a tridiagonal system which can be solved using either 
alternating direction implicit (ADI) or successive over relaxation (SOR) 
methods. The reduction of the equations to tridiagonal form has been 
developed in Appendix B. 

To implement an ADI solution, a Thomas algorithm [52] is first 
applied along each line of constant jk for the x direction sweep. This 
gives the three components of vorticity at the centers of the x=constant 
cell faces. Explicit in this sweep are the three components of 

vorticity occurring on the y=constant and the z=constant faces. For the 
first iteration through x direction sweeps, the values of the vorticity 
components are at time level n. The Thomas algorithm is then applied in 
the y direction and z direction, resulting in the updated components of 
vorticity on the y=constant and z=constant faces respectively. In 
general, it is found that sufficient refinement of the solution is 
achieved by cycling through each sweep direction three times. 
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Recall that the velocity coefficients of the advection terms are 
lagged by one time level during these cycles. To achieve second order 
accuracy in time it is necessary to update these velocities to the 
current time level of the vorticity, and then recompute the vorticity 
using these updated velocities. It is also necessary to recompute the 
right hand side of Eq. (3.49) after recycling through the velocity 
equations and before recomputing the vorticity. This requirement 
results from the fact that the n+1 time level of the exponential 
transformation appears in Eq. (3.49). Repeated application of the above 
procedure yields in the vorticity and velocity fields at later times. 

Alternately, the implicit system may be solved using a SOR type 
iteration scheme instead of applying the Thomas algorithm along lines. 
The advantage of this approach is that when only a few cycles are 
required, the SOR iteration will be faster computationally. In 
addition, a residual can be identified and used as a criterion for 
convergence. The residual is defined as the difference between the left 
side and right side of each individual equation of the tridiagonal 
system. The other aspects of the solution procedure are identical with 
the method employing the Thomas algorithm. In all the subsequent 

computations, the time steps were kept small enough to ensure that the 
CFL stability criterion [53] was not violated. 


Helmholtz Projection 


In general, the vorticity resulting from the finite difference 
solution to the vorticity transport equations does not satisfy the 
requirement that the vorticity vector be divergence free. This is due 
to the fact that the divergence free condition is a vector identity and 


42 


is not part of the Navier-Stokes equations. Hence, the divergence free 
condition is not required in the derivation of the vorticity transport 
equation. Thus, an initially divergence free field may drift from the 
requirement due to roundoff error or discretization error. 

A previously discussed requirement for the numerical solution of 
the velocity field was the compatibility condition 7*u=0. If the vector 
resulting from the time advance of the transport equation does not 
satisfy this condition, it must be projected onto the vector space of 
divergence free vorticity. A well known property of any vector field is 
the fact that it can be decomposed into an irrotational field and a 
divergence free field according to the Helmholtz projection: 


<*> = ?X + 7 x c (3.58) 

irrotational divergence free 

Here, x represents a scalar function and c , a vector function. To 
extract the divergence free part of the vector, <j, take the divergence 
of Eq. (3.58) to get: 

7»« 3 7 2 x ' (3.59) 

Equation (3.61) can be used to solve for x . subject to the boundary 
condition ~ * 0 on B. This boundary condition is utilized to ensure 
that the component of vorticity normal to the boundary is not altered. 
The divergence free component, 7 x c, of the original vector, Z, is 
given by 

7 x c 3 u - 7 X (3.60) 
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An alternate approach to enforcing the solenoidal condition on the 
vorticity vector can also be derived. Equation (3.59) can be reduced to 
the following set of first order differential equations: 


H*$ ?♦&■*•“ 

(3.61) 

p -f 

(3.62) 

a - dX 

Q " ay 

(3.63) 

r - IS 
dz 

(3.64) 

where ~ = 0 on the boundary. The unknowns in 
on 

this system are the 

scalar quantities p, q, r and x. The divergence 

free components of the 

original vector are then given by: 


(7 X C)j * (Jj • p 

(3.65) 

(V x c) 2 * - q 

(3.66) 

(7 x c) 3 * - r. 

(3.67) 

The di scretization of equations (3.61) to 

(3.64) is described by 

Rose [54]. The finite difference forms are: 


5 p + 5 q + 6 r = 7*u 
x y^ z 

(3.68) 

“x p = V 

(3.69) 

M« y q * s y x 

(3.70) 

v-r - 

(3.71) 
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(3.72) 


V ‘ 5 x p = V ' V 
“x* ' TT 5 x p * “z* ■ T~ V (3J3) 


The averaging and difference operators are identical to those defined 

previously for the discretization of the transport equation. 

These equations have been solved using the Kaczmarz [48] iteration 

scheme which was discussed in reference to the solution of the velocity 

field. Recall that the Kaczmarz scheme requires the evaluation of the 
-T - -T -1 

expression A. (A. A.) R_. for each of the discretized equations, the 

subscript i denoting a specific equation. The expressions for A^ and 
R. are given below for Eqs. (3.66) zo (3.7 i) as tqs. (3.74 a.b) to (3.80 
a,b) respectively. 


a , (i. li L i L ii) 
*1 'Ax* Ax* Ay* Ay* AZ AZ 


(3.74a) 


R 1 “ Ax ^ P i+V2 - P i- V 2 ^ + Ay ^ j+ V 2 ” q j- V 2 


1 — 

+ — (r, 1 - r 1 . ) - (7»aj) . . 

Az k+ V2 k-V2 i.J.k 


(3.74b) 


A . (I 1 li L.) 
2 2* 2’ Ax’ ax 


(3.75a) 


' 2 3 T {p i + v ? +p i-W 


TZ ( X i+ V 2 ‘ X i- V2 


(3.75b) 


A * (i i ZL L) 
3 2* 2* Ay’ Ay 


(3.76a) 
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(3.76b) 


R 3 ’ ? (q j + V 2 + q j- V 2 1 - W <X J> l/ 2 ' X j- V 2 ’ 

, M i I I!, 

4 2* 2’ Az- Az' 

R 4 *1 tr k*V W 2 1 "k ( W 2 ‘ Vl/2 1 

.,11 _1 1 -Ax Ax Ay -Ay. 

5 " 2’ 2’ * 2* ’ V 3 ’ 8 ’ 8“* ~ 

r 5 ■ 1 ( x i + V 2 * x i- V 2 ’ * 7 ( V 2 + x j- v 2 > 

AX 2 2 

- -r ( p, + v 2 • p,-- v 2 1 + T- 1 v v 2 - v v 2 > 

. 3 4 l 1 1 AX AX AZ AZ. 

6 2 ’ 2 ’ “ 2 ’ ‘ 2 ’ ~ 3 ’ 8 * 8 * ~ 8 ~ 

r 6 ? <*i+V X <-V 2 ’ • ? ( \.V 2 * W 2 ’ 

4x 2 2 

" 5~ (p i+ 1/ 2 ' p l- V 2 1 * T (r k+ V 2 " r k- V 2 ’ 


Each of the above equations is relaxed independently over all 
and boundary cells using the Kaczmarz iteration scheme repeated 






T -1 


(A. A.) 


R. (k) 


(3.77a) 

(3.77b) 

(3.78a) 

(3.78b) 

(3.79a) 

(3.79b) 

interior 

below. 

(3.80) 
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The vector it is composed of combinations of the variables p,q,r, and 
X depending on the equation being relaxed. For instance, when i=l, 

which corresponds to Eq. (3.68) the vector n, is given by 

*1 ‘ (p i+v 2> P >-V 2 . W 2 . p j-v 2 , VV 2 . r k-v 2 .> (3 - 81) 

The boundary conditions are handled in a similar manner. It is 

required that the normal component of vorticity on a boundary remain 

unchanged after the Helmholtz projection. Thus, on constant x boundary 

faces ~ = 0, on constant y boundary faces, -r 8 0, and on constant z 
gX oy 

^ V 

boundary face, ~ = 0. Hence, the boundary conditions are 


p = 0 

on 

i-1, and 

imax 

(3.82) 

xs 

II 

o 

on 

j»l, and 

jmax 

(3.83) 

r = 0 

on 

k=l , amd 

kmax 

(3.84) 


These values are reset after each iteration through the Kaczmarz update, 
since in general, the boundary values of p,q and r change after the 
relaxation of Eqs. (3.68) to (3.73). 

An alternate method, mentioned initially, is to solve the Poisson 
equation for x» and then compute Vx. This can be done using an SOR 
iteration scheme. The divergence of vorticity is given at the center of 
a cell, so that x should be computed at cell centers while *?x should be 
computed at cell faces. To facilitate coding, it is desirable to 
express the Poisson equation in terms of the orthogonal, curvilinear 
coordinates 

x = x (£) 
y = y (n) 
z = z (c) 



The transformed Poisson equation takes the form 





XZ 





(3.88) 


where 

0<C<1, 0<C<1 


The computational plane consists of rectangular cells with uniform 
spacing in each coordinate direction. This allows one to use simple 
unweighted finite difference expressions. Equation (3.85) is 
discretized using central differences for both first and second order 
derivatives. The discretized form can be solved using an SOR scheme 
for the variable x • It is then necessary to compute Vx as follows. 


?x 



+ n j + c -22 k 

y dr J ^z dr 


(3.86) 


In the discretized form, the first order derivatives are represented by 
central differences. On boundaries, dx/3n = 0 , so Eq. (3.86) is 
computed for interior cells only. If the divergence free components are 
represented by then 

«'l ■ «1 “ 3x/dx * (3.87) 
w' 2 3 u 2 “ (3.87b) 
u)' 3 * u 3 - dx/az (3.87c) 


Pressure Equations 

One advantage of the vortici ty-veloci ty formulation is that the 
pressure does not appear explicitly in the equations of motion. 
Therefore, the velocity and vortici ty are obtained independent from the 


48 


pressure field. In Appendix C, it is shown that pressure satisfies a 
Poisson equation, whose right hand or function side is an expression 
containing the velocity gradients. The numerical solution of this, 
equation is analagous to the solution of the Poisson equation, £q. 
(3.59). 

Neumann conditions, resulting from the momentum equations, were 
specified on all boundaries. A special requirement of these boundary 
conditions, due to an application of Green's theorem, is discussed 
below. 

The integral form of the Poisson equation for pressure can be 
developed from the momentum equations and written as: 

j T^P dY = j Sdv (3.88) 

where the source term S is given as 

_ r5v au aw au aw av au av au aw av aw-i , , 

b = ~^ p '■ax ay ax az ay az ax ay “ ax az ay az J 1 * syj 

Through an application of Green's theorem, a relationship between the 
source term and the boundary flux is given as: 

/ SdV = / aa (3.90) 

Here, n is taken as positive when directed outward from the boundary 
a. In general, the finite difference equivalent of Eq. (3.90) will not 
be satisfied identically. As a result, the numerical solution of the 
Poisson equation for pressure need not converge. 

A convergence requirement can be developed for the Poisson equation 
by imposing an error condition, E, defined as 

E = / SdV - / ~ aa (3.91) 
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The strategy is to distribute this error over the boundary flux terms 
dP/dn, which represent the boundary conditions. In this sense, the 

finite difference analog of Eq. (3.90) is satisfied. 
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SPECIFICATION OF THE PROBLEM - 
80UN0ARY CONDITIONS 


Mathematically, the Navier-Stokes equations are a set of three 
elliptic, second order partial differential equations. This means that 
either a Dirichlet or Neumann condition must be specified on all closed 
boundaries. It is these boundary conditions that distinguish all the 
different flow patterns occurring in nature. Therefore, it is extremely 
important that these conditions be chosen properly. 

The purpose of the present study was to model the evolution of a 
class of vortices similar to those shed from the wing tips of 
aircraft. Although these vortices generally occur in pairs, only a 
single vortex was considered. Observations reveal that these vortices 
are roughly axisymmetric with an appreciable axial velocity (either a 
defect or an excess velocity relative to the free stream, depending on 
the wing loading [55]). In. addition, far downstream of the wing the 
flow outside of the core is nearly irrotational . 

Experimental measurements of a trailing wing tip vortex [42] and 
vortices produced in a tube and vane apparatus C27 ] revealed that the 
circumferential velocity profile is well represented by the two- 
dimensional 8urgers vortex. The axial velocity appears - to decay 
exponentially in the radial direction (from the vortex centerline), 
reaching a constant value at large radius. These observations can be 
represented by the following dimensional form for the swirl and axial 
velocity profiles. 


2 

-ar 

v = ( K/r) ( 1-e 2v ) (4.1) 

e 

2 

-ar 

u = U + U e 2v (4.2) 

« o 

Here, K is a constant, which is proportional to the circulation, "a" is 
an adjustable constant associated with the vortex core diameter and v is 
the kinematic viscosity. is the free stream axial velocity and U Q is 
an axial velocity excess (or deficit) occurring at the vortex center- 
line. Profiles similar to these were employed in numerical studies of 

vortex breakdown by Kopecky and Torrance [31] and 8enay [35]. 

In the above equations, length was nondimensiona 1 ized by r' , where 

r' = /2Va (4.3) 

and velocities were nondimensiona! ized by the axial velocity, U*, 
occurring at the radius of maximum swirl velocity, r*. Formally, the 
quantity r'is the distance in which the axial vorticity e-folds once. 
In fact, the radius of maximum swirl velocity- is nearly equal to r*. 
Through an .iterative process one obtains the relation 

r* = 1.12 r' (4.4) 


The nondimensiona! forms of Eqs. (4.1) and (4.2) become 


1 1 ,, -r , 

v 9 s rrrto? (1 - e 1 


(4.5) 


u = 


1 + 6e r 


1 + 0.2856 


(4.6) 


★ * _ 

where the Rossby number is given as Ro = U /r Q and 6 = U /U . Here. 

0 oa ’ 

Q is defined as the solid body rotation rate obtained from 

2 * lim (V/r) = aK/2v. (4.7) 

r-o 

The axial component of vorticity, for the velocity distribution given by 
Eqs. (4.1) and (4.2), is given as 




X 



(4.8) 


By defining Q as the reference vorticity, the above equation can be 
written in nondimensional formas 


-r 


“x * 2 e 


The circumferential vorticity component is given as 

.2 


ar .. 

"9 = ~ V 


which in nondimensional form becomes: 

2.246 


-ar 

77 


U 9 


Ro re 


-r 


1.0 + 0.2856 


(4.9) 


(4.10) 


(4.11) 


These expressions represent the form of the vortex employed in the 
present study from which the initial and boundary conditions were 
derived. The initial and boundary conditions at inflow are discussed 
first, followed by the outflow conditions and conditions at large 
radius. 

The specification of the velocity vector, Eq. (3.3), or its 

gradient on all closed boundaries is required to solve the primitive 
variable form of the Navier-Stokes equations, Eq. (3.2). The numerical 
solution of the vorticity transport equation, Eq. (3.5), requires the 
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specification of three components of vorticity on the boundaries. The 
numerical solution of the “velocity equations", Eqs. (3.1) and (3.4), 
requires the specification of a single component of the velocity vector 
(or gradient) on the boundaries. The general procedure is to specify 
the primitive variables on a boundary and translate this specification 
to a finite difference form consistent with the above requirements of 
the algorithm. 

Although qua si -cylindrical vortices are best described in a 
cylindrical coordinate system, the algorithm, as constructed at this 
point, was written in terms of Cartesian coordinates. Therefore-, 
references to both coordinate systems are required in order to interpret 
the influence of the boundary conditions on the solution. In the 
following discussion it is assumed that the initial vortical 
distribution was cylindrical (i.e., no variations in the axial direction 
and no radial velocities). The vortex was aligned along the x axis and 
the rotation was such that the axial component of vorticity was 
positive. 

At inflow, the specified primitive variables were: the axial 

velocity component, u, the y derivative of the w, and the z derivative 
of the v components of velocity, (i. e. u, — , -^y) . These conditions 
must be interpreted with respect to the algorithm. The axial velocity 
component, u, was held fixed at inflow when solving the "velocity 
equations." This allowed for direct specification of wake-like or jet- 
like profiles. The v and w components of velocity were specified to the 
extent that they took on the values determined by the component of 
vorticity normal to the inflow boundary. This led to the specification 
of dw/by and 3v/dz as the additional primitive variable boundary 
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conditions, since these derivatives can be combined to give the axial 
component of vorticity. Note further that the axial component of 
vorticity can be oriented in the same direction in both cylindrical and 
Cartesian coordinate systems. Therefore, the specification of a> in the 
Cartesian coordinate system (of the algorithm) can be interpreted as 
specifying the vorticity distribution of Burgers’ vortex, since dw/3y 
and 3v/3z were obtained analytically from Eq. (4.5). 

The "vorticity solver" required the specification of two additional 
boundary conditions. These were chosen as u) y and oj^. 

In nondimensional form, u» and u> can be written: 

y z 


“y * U2 Ro 'a? ' a? 1 


(4.12) 


«z = U12 Ro " $ 


(4.13) 


Since u is an analytic boundary condition at inflow, the derivatives 
au/dz and 3u/3y were known. It was necessary to calculate the deriv- 
atives aw/ax and av/ax numerically since they could not be derived from 
the analytic boundary conditions. First order forward differences were 
used to compute these derivatives. The resulting boundary conditions 
were assumed to be at time level n (present time). 

A different strategy had to be used to specify the outflow boundary 
conditions since the solution here was unknown and highly dependent on 
the flow upstream. With respect to the "velocity solver", the condition 


av/ay + aw/az = constant 


(4.14) 


was chosen. This is a statement regarding the divergence of the 
velocity in the plane perpendicular to the vortex axis. The constant 
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equals zero for a steady flow in the limit as the Rossby number 
approaches zero [40]. 

Flux conditions were chosen for the vorticity boundary conditions 
at outflow. Here, the effects of viscosity were neglected and the 

Dirichlet type boundary conditions were obtained assuming a time advance 
of 

(4.15) 

The velocities and vorticities on the right hand side, of Eq. (4.15) were 
taken at time level n. The time derivative was discretized using first 
order forward differences. The resulting discretized equation was then 
solved for the vorticity vector at time level n+1. 

At the large radius boundaries given by planes of constant j and k, 
the axial component of velocity was specified. This was done so that 
the effects of an external pressure gradient, analagous to the 
experimental investigations, could be modeled. 

At these boundaries the three components of vorticity were 
specified using the Cartesian coordinate equivalents of Eqs. (4.9) and 
(4.11). For a uniform inflow profile, Eq. (4.11) shows the and 

<*> 2 components of vorticity are zero. By evaluating Eq. (4.11) at large 
radius, one can see that the axial component approaches zero asymptot- 
ically. Therefore, the radial boundaries were placed at a radius which 
was large enough to ensure that vorticity did not, through convection or 
diffusion, contaminate the boundary conditions. 

To summarize, the specification of the boundary conditions for both 
the "velocity solver" and the "vorticity solver" are represented in 
Table 4.1. 
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Table 4.1 


Summary of Velocity and Vorticity Boundary Conditions 


Velocity Solver 



Surface 

Speci fica tion 


(a) inflow 

u=gi ven 


(b) outflow 

3v/3y + dw/az = 0 


(c) radial boundaries 

u=gi ven. 


Vorticity Solver 



Surface 

Speci fica tion 


(a) inflow 

- r 2 

“x * 2 e 

co 3 1.12 Ro (3u/3z 
y 

- aw/ax) 


<o z 3 1.12 Ro (5v/3x 

- 3u/3y) 

(b) outflow 

Deo - - 

— = U *V<o 
Dt 


(c) radial boundaries 

io 3 given 




The discretized form of the governing equations were solved over a 
48x28x28 grid (47 cells in the x-direction, 27 cells in the y-direction, 
27 cells in the z-direction) on a Cyber 205, Course grid (52x20x20) 
solutions were computed on a CDC 830 at Old Dominion University to help 
identify relevant parameter ranges. Since the difference scheme is 
compact, grid clustering is easily performed. Grid points are usually 
clustered in regions where large gradients occur. In studies of vortex 
breakdown, this includes the region immediately upstream of the 
breakdown and the vortex core region. Immediately upstream of 
breakdown, large gradients of the axial velocity occur as the fluid 
approaches the stagnation point. Large gradients of the circumferential 
velocity are present within the core. 

Refinement of the mesh in the axial direction was accomplished 
using the transformation: 


x 


h (o 2 -l) [(-—) - 1] 
x 3-1 

x-l 

(a*n ♦ i] 

a-1 


(4.16) 


where 3 is a stretching parameter and h x is the length of the domain in 
the x direction. More points are clustered near x=0 as a+1, l<a<®. 
The coordinate x varies uniformly from 0 (corresponding to the x=0 grid 
point) to 1 (corresponding to the x=h x grid point). 

To cluster points near the vortex centerline (y and z coordinates) 
the following transforma tion was used: 


y=y {l 


^ sinhCe(y-B) ] , 


sinh ( e8 ) 


(4.17) 



where 


. 1 + (e -1) (y /h ) 

[ - M-1. 0<=< 

1 + (e -l)(y /h ) 

c y 


and e is a stretching parameter. Increasing e clusters more points near 
y=y c . If e=0, a uniformly spaced grid results. The domain length is 
given by h^. The coordinate y varies uniformly from 0 (corresponding to 
the y=0.0 gridpoint) to 1 (corresponding to the y=h^ gridpoint). 

The above transformations can be found in Anderson, Tannehill and 
Pletcher [52], The values of the grid parameters for the solutions 
generated in this work were given as: 

h^ * io.O, a - i - 15 x direction grid 
y * 5.0, h = 10.0, s = 4.5 y direction grid 
z c = 5.0, h z = 10.0, e = 4.5 z direction grid 

These resulted in the following minimum cell lengths: 

Ax * 0.1303 
Ay » 0.1778 
Az * 0.1778. 


The distribution of the cells within the domain and the orientation 
of a typical vortex is shown in Fig. 4.1. 
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Fig. 4.1 Schematic representation of a typical vortex within the 
computational domain. 


RESULTS AND DISCUSSION 


Results of two test cases are presented. In the first case, 
(Ro=0.625, Re=225), the vortex was imbedded in a uniform external flow 
with no external pressure gradient. The flow produced by this Rossby, 
Reynolds number case, along with numerous other cases, predicted an 
axi symmetric breakdown which occurred at the inflow boundary. That 
result was similar to the results obtained by previous investigators 
[31,33,35,36]. The questions which arise from breakdown near inflow 

have been discussed earlier in this work. In order to alleviate this 
problem, a different type flow was computed, (Ro=0.8, Re=225), in which 
the vortex was imbedded in a decelerating free stream. This, in 
essence, modeled the effects of an adverse pressure gradient on the 
streamwise development of the vortex. The resulting breakdown occurred 
away from the inflow boundary. In addition, a multiple celled breakdown 
region was observed, in accord with experimental observations [27]. In 
both cases, the parameter, 6, defined in £q. (4.6), equaled 0.0. Thus, 
a uniform axial velocity profile was specified at inflow. In the 
following discussion, detailed results of the above two test cases are 
displayed in the form of line and contour plots. Several other test 
cases that were run without graphic output are discussed. Finally, the 
results of tests performed to ascertain the effects of domain length and 
grid size on the temporal evolution of the solution are discussed. 
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Case i - Uniform Free Stream Axial Velocity 


,3ased on previous numerical studies, breakdown of trailing wing tip 
vortices occurred at Rossby numbers of approximately 0.6 or less when 
the Reynolds number was greater than 100. This behavior is shown 
clearly in Fig. 2.1. The initial numerical simulation in the present 
study was completed in an attempt to verify these results and to 
determine the effects of asymmetries on the solution. The test was 
performed at a Rossby number of 0.625, which was chosen because it was 
near the delimiting line for vortex breakdown displayed in Fig. 2.1. 
Flows computed by previous investigators at Rossby numbers considerably 
below the delimiting line in Fig. 2.1 appear to become distorted and 
non-physical near the inflow boundary. A Reynolds number of 225 was 
chosen to minimize the apparent damping effects of viscosity at very low 
Reynolds numbers. The results of this simulation are displayed in Figs. 
5.1 to 5.15 for the time level t=126.8. Unless otherwise noted, the 
contour plots are in the x-y plane along the centerline of the vortex. 
Solid contour lines denote positive values (or zero) and dashed lines 
denote contours with negative values. In all plots, the contour levels 
are evenly spaced. Plots of particle traces, vortex lines and velocity 
vectors are also displayed. They were obtained using PL0T3D, a three- 
dimensional color graphics program implemented on an Iris color graphics 
workstation. These plots are projections of three-dimensional vector 
fields onto a two-dimensional surface. 

Velocity vectors, projected onto the mid plane (z=5), are displayed 
in Fig. 5.1. Observe that the internal structure of the breakdown 
region is seen to consist of a single cell, nearly symmetric about the 
vortex centerline. Fluid is entrained through the top half of the rear 
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Fig. 5.1 Projected velocity vectors over the Interior of the breakdown 
region, (x-y plane) 


of the breakdown region. The fluid appears to exit the bubble fron the 
lower half. In addition, two stagnation points along the axis can be 
di stinguished. 

Figure 5.2 . represents particle traces. Nine white colored traces 
were started at the inflow plane. These traces were all started from a 
radial position within the rotational portion of the core. In addition, 
a trace was started within the breakdown region itself, and is 
represented by the red line. A third trace was started at the inflow 
plane but outside the core in the irrotational region of the flow. That 
trace is blue. Particle traces satisfy the equations dx/6t=u(x, t) . If 
a particle passes through the point (x,y,z) at time t=0 the solution is 
of the form x=x(x,y,z,t) which traces out the pathline as t increases. 
The PL0T3D graphics package is limited to instantaneous particle traces, 
i.e., the velocity components must be time independent. Therefore, in 
this sense, the traces can be considered as streamlines, pathlines or 
streaklines because in a steady flow they all coincide. The tangents to 
these traces are everywhere parallel to the velocity vector. 

The white traces shown in Fig. 5.2 can be seen to approach the 
breakdown region and diverge - never entering the cell itself. These 
traces also reveal that the diameter of the vortex core has increased 
behind the breakdown region. The red trace, released from within the 
bubble, is seen to spiral about within a single cell. It eventually 
exits the region from behind, closer to the axis radius than the white 
traces. The blue trace spirals around the outside of the vortex core in 
the irrotational region and is essentially unaffected by the breakdown. 

Contours of the axial velocity are displayed in Fig. 5.3. This 
plot clearly reveals a breakdown region which occurs in close proximity 
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Fig. 5.2 Instantaneous streamline pattern for flow with vortex 
breakdown. 


to the inflow boundary. The recirculation region is bounded by the area 
within the innermost solid contour line (the 0.0 level contour). The 
breakdown bubble appears symmetric about the vortex centerline. 

Figure 5.4 represents the w (axial) component of vorticity. For 

this figure and all succeeding vorticity contour figures (including w , 
and w contours), the vorticity has been scaled by the Rossby number. 

For example, u> x = 1.12 Ro(dw/3y-3v/3z) . An interesting feature of this 
flow is the intensification of the axial component of vorticity 
occurring just aft of the breakdown region. Here, the axial vorticity 
has increased by approximately 252 over the maximum value at inflow-. 
This is due to vortex stretching which results from a rapid acceleration 
of the axial velocity component. Within the breakdown region itself, 

the axial component of vorticity is small. This indicates that within 
the breakdown region the radial gradients of the circumferential 
velocities are small. Downstream of the breakdown region the contour 

lines become parallel, revealing a return to a quasi-cyl indrical flow. 
In the absence of breakdown, the vorticity contour lines over the entire 
region would appear nearly parallel. 

Contours of co y and w 2 vorticity are displayed in Figs. 5.5 and 5.6, 
respectively. The 0.0 level contours are not displayed, because outside 
the core of the vortex, large regions exist where the vorticity field is 
near zero. Therefore, plotting zero level contours results in a large 
number of undesirable and confusing lines in the far field. The plot 
represented by Fig. 5.5 is in the x-z plane since the component of 
vorticity in the x-y plane is nearly zero. The w and components are 
due entirely to perturbations of the base flow since the initial (t=0) 
distribution possessed only an axial component of vorticity. 
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Fig. 5.4 Contours of constant axial vortlcity, Contour levels 
range from 0.25 to 2.5 in intervals of 0.25. 
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Fig. 5.6 Contours of constant w 2 vortlclty. Contour levels range from 
-1.25 to 1.25 In intervals of 0.25. 
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Vortex lines are shown in Fig. 5.7. These are lines whose tangent 
is everywhere parallel to the vorticity vector. The lines are three- 
dimensional and were obtained using the Iris color graphics work-station 
and PL0T30. As was the case with the particle traces, the vortex lines 
must be considered to be time independent. The radial locations of the 
lines at the inflow plane correspond to the radial locations of the 
white particle traces in Fig. 5.2. Note that the magnitude of the 
vorticity cannot be inferred from these lines. At inflow, the lines are 
oriented in the x direction, revealing that only the axial component of 
vorticity exists. Upon entering the breakdown region the lines become 
oriented nearly perpendicular to the vortex axis. This signifies a 
transfer of vorticity to the w , and components, and is controlled 
by the vortex stretching and bending terms in the vorticity transport 
equation. Aft of the breakdown region, the vortex lines are oriented 
primarily in the axial direction, indicating a transfer of vorticity 
back to the axial component. 

1 -2 

Figure 5.8 is a line plot of the integral J j u> dV as a function 
of axial location, where u> /2 is defined as the enstrophy. The 

integration was carried out over control volumes defined by 0<y<10, 
0<z<10, Ax. The resulting numerical values were then divided by the 
volume over which the integration was ' performed. When the enstrophy is 
integrated over a volume of fluid it is an appropriate measure of the 
total amount of vorticity within the fluid [40]. The enstrophy is 
maximum within the breakdown region at the axial location x*3. The 
enstrophy decreases and remains nearly constant downstream of the 


bubble. 
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Figure 5.9 is a line plot of the volume integral of the material 

derivative of enstrophy as a function of axial location. The integral 
d 1 -2 

is given as ^ / j w dv which provides a measure of the rate of change 
of the total amount of vorticity in a specific volume of fluid, V. 
Here, V is enclosed by a surface 3 moving with the fluid. The plot 
displays axial variation of total enstrophy, along with the rate of 
change provided by stretching, viscosity and the flux of enstrophy 
across the boundaries of the control volume. The meaning of these terms 
can be explained by expanding the time rate of change of enstrophy 
integral in terms of volume integrals containing Eulerian derivatives-. 
The result, derived in Appendix 0, can be written as: 


rl 1 bu . dw . . 

3t ! 2 * I "i“j ST dv ■ ‘-STT’ dv * 2 


3(cj.m. ) 

t 11 

i — 


dS (5.1) 


By examining the right hand side of Eq. (5.1) it can be seen that the 
total amount of vorticity in a material volume can change as a result of 
vortex stretching and viscous effects. The first term on the right, the 
stretching term, is positive if the fluid element is extended in the 
direction of the local vortex lines. The second term reveals that the 
effect of viscosity, neglecting diffusive transport across the 
boundaries, is to decrease the total enstrophy of the fluid. In 
general, it is possible for the entire right hand side of Eq. (5.1) to 
be positive, leading to an increase in the enstrophy, or total amount of 
vorticity, in the fluid. 

The volume integrals on the right-hand side of Eq. (5.1), which are 
plotted in Fig. 5.9, were evaluated using trapezoidal rule integra- 
tion. The integrations were performed over volumes defined by 
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Fig. 5.3 Variation of Integrated enstrophy with axial location. 
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Fig. 5.9 Variation of integrated material derivative of enstrophy 
with axial location. 
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0<y<10, CKz<10, Ax. The plotted values represent volume averages, i.e., 
the numerical values resulting from the trapezoidal rule integration 
were divided by the volume over which the integration was performed. 

From a lagrangian point of view the rate of change of enstrophy of a 

material volume is. due to both temporal changes and changes due to 

spatial movement of the volume. In a steady flow, the temporal changes 
are non-existent. Since the breakdown is an unsteady phenomena, 

temporal changes may be significant although it is unlikely that they 
are dominant. Whenever the rate of change of enstrophy is negative, the 
total amount of vorticity contained in a material volume passing through 
that location is decreasing. By examining Fig. 5.9 it is apparent that 
the distribution of enstrophy within the fluid is controlled, for the 
most part, by vortex stretching. Viscous effects appear to be important 
only within the breakdown region. The diffusion of enstrophy into, or 
out of, a material volume is negligible. The maximum and minimum values 
of the total rate of change of enstrophy occur at axial locations 
corresponding to the stagnation points. These points approximately 
define the front and rear of the bubble. 

* 2 

A contour of the pressure field, non-dimensional i zed by pU , is 
shown in Fig. 5.10. Beginning at the inflow boundary and near the 
centerline, the fluid encounters an adverse pressure gradient with a 
corresponding decrease in the axial velocity. Within the breakdown 
region itself, the pressure remains nearly constant. The fluid is 
accelerated beyond the breakdown region, and this is manifest in the 
form of a weak favorable pressure gradient. Oownstream, the pressure 
contours are nearly parallel. Since the vortex was imbedded in a free 
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stream with a constant axial velocity, the pressure along the radial 
boundaries is nearly constant. 

The relationship between the centerline axial velocity and the 
centerline pressure is shown in Fig. 5.11. Except for a short distance 
within the breakdown region, increasing axial velocities correspond to 
favorable pressure gradients and vice versa. 

Figure 5.12 is a contour plot of the dimensionless total pressure, 

? 2 

q /2 + p. In a steady inviscid flow, q /2 + p remains constant along a 

streamline. Therefore, the above contours can be thought of as 

approximating streamlines. Along a streamline, wherever the dynamic 
2 

pressure (q /2) is high, the static pressure (p) is low and vice versa. 

The rate of change of energy is plotted in Fig. 5.13. The rate of 
change of total energy is given as the sum of the rates of change of 
internal and kinetic energy. The rate of change of kinetic energy is 
given as the sum of the rates of change due to both pressure and viscous 
forces. For an isothermal, incompressible fluid in a lagrangian 
framework, the rate. of change of total energy of a specific material 
volume is given as: 


d 

Tt 


I (E 


+ j u 2 )dV = 


-/ — dV ♦ v/u. 

J p ax, J i 


aV 


dx • dx . 

J J 


dY 


+ dv/e . .e. . dV (5.2) 

' J ■ J 


The interpretation of the above material derivative is analagous to the 
interpretation given earlier for the enstrophy equation. The rate of 
change of internal energy is due solely to the dissipation of mechanical 
energy through viscous effects. Here, the shape of the fluid element is 
distorted. This is an irreversible loss of energy, manifested in the 
form of heat. This is accounted for the third term on the right hand 
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VELOCITY AND PRESSURE 



Fig. 5.10 Isobar plots: contour levels varying from -0.9 to 0.1 in 
intervals of 0.1. 
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Fig. 5.11 Variation of axial velocity and pressure along the vortex 
centerline. 
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Fig. 5.12 


Contours of total pressure: contour levels of -0 
in intervals of 0.1. 


.5 to 0.7 
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side of Eq. (5.2). The effect of the first and second terms (pressure 
and viscous stresses, respectively) is to change the kinetic energy of a 
fluid element. The pressure forces are reversible changes in energy. 

The kinetic energy per unit volume is plotted as a function of 
axial distance in .Fig. 5.14. Following a slight decrease near inflow, 
kinetic energy remained constant in the axial direction. 

Figures 5.15 a-d represent velocity profiles at four different 
axial locations. Axial (x), transverse (y) and spanwise (z) velocities 
are plotted as a function of the transverse coordinate (with a 
translation of the origin to the vortex centerline; y=5.0, z=5.0). 
Since the spanwise location of the data points was along the vortex 
centerline, the velocity components in a Cartesian system can be 
transformed easily to the corresponding components in a cylindrical 
system. The axial location of the profiles in Fig. 5.15a is slightly 
upstream of breakdown (x=0.27). Figure 5.15b represents profiles from 
within the breakdown region (x=2.52). Figure 5.15c represents profiles 
near the aft stagnation point (x=3.85). The profiles plotted in Fig. 
5.15d are at an axial location downstream of the breakdown region 
(x=9.52). Upstream of breakdown, the spanwise (swirl) velocity 
represents the two-dimensional Burgers vortex to a good approximation. 
The flow is no longer quasi-cyl indrical , since non-zero transverse 
(radial) velocities appear near the centerline. These velocities are 
due to the divergence of the flowfield away from the axis as the 
stagnation point is approached. The axial velocity shows approximately 
a 20% deficit near the axis. At an axial location within the breakdown 
region, decay of the spanwise (swirl) velocity at large radii (outside 
the recirculation zone) is inversely proportional to the radius. Within 
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Fig. 5.13 Contribution of various energy transport terms to the 
integrated rate of change of energy. 



Fig. 5.14 Variation of integrated kinetic energy with axial location. 
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Fig. 5.15 Radial variation of axial, swirl and radial 
a,b,c,d. components at different axial locations 
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the recirculation zone, a short annular region is observed in which the 
spanwise (swirl) velocity is constant. The transverse (radial) 
velocities are near zero for all transverse (radial) locations. 

In Fig. 5.15c, the transverse (radial) velocities have become 

significant, and are nearly the same magnitude as the transverse 
(radial) velocities in Fig. 5.15a, but are of opposite sign. This 
signifies that the axial location is near the aft portion of the bubble, 
since the streamlines are converging. 

It is revealed in Fig. 5.15d that the flow downstream of the 

breakdown is quasi-cyl indrical , with a small axial velocity deficit. It 
is apparent that the vortex core diameter, defined by the radius of 
maximum swirl velocity, is greater downstream of the breakdown region 
than upstream. 

The minimum axial velocity varies with time as shown (on a semi-log 
plot) in Fig. 5.16. Note that the axial velocity decays exponentially 

to approximately 30 percent of its original magnitude (over a time 

period of 60 units). An exponential decay of the velocity field is 
indicative of a viscous time scale. By dimensional analyses, the 
viscous time scale for the flow is given as t = 1.12 Re, which for this 
case, equals approximately 100 time units. 

Case 2 - Deceleration of the Free Stream 
Axial Velocity 

The breakdown just described was very similar to breakdowns 
computed in several previous investigations. In an attempt to alleviate 
the problem of the non-physical results and breakdown occurring near the 
inflow boundary, an additional vortex flow was computed. Here, the 
vortex was imbedded in a decelerating free stream. The purpose of this 
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case was to force the local Rossby number to decrease as the flow 
evolved in the axial direction. The Rossby number at inflow was 0.8 and 
the Reynolds number was 225. The free stream axial velocity was 

linearly decelerated from 1.0 to 0.55 over the range x=1.43 to x=16.0, 
thus du/dx=0.03. The results are plotted in Figs. 5.17 to 5.32. The 

data in Figs. 5.18 to 5.32 is at time level t=81.28. Five different 

time levels are represented in Fig. 5.17. 

Velocity vectors, representing time levels t=81.28, 85.27, 87.45, 

89.63 and 91.82 are displayed in Fig. 5.17 a-e, respecti vely . The 
general appearance of the bubble is one of asymmetry. At time t=81.28 
the internal structure of the breakdown contains two major cells, or 
vortex rings, rotating in opposite directions about their respective 

axis. The aft vortex ring is inclined to the x-axis. Fluid enters the 

bubble from near the downstream end, through the side of the bubble, and 
exits the bubble at the same axial location but on the opposite side. 
The inclination of the aft vortex ring also appears to be related to the 
exchange of fluid in the bubble. The most forward section of the ring 
corresponds to the location of fluid influx, whereas the aft section of 
the ring corresponds to the location where fluid is emptied. In 

addition, the velocities are considerably greater in the aft portion of 
the bubble than in the forward portion. The length to diameter ratio of 
the bubble is approximately 1.75. The maximum diameter of the bubble 
occurs approximately 0.7L units downstream from the front stagnation 

point (where L is the length of the bubble). 

The velocity vectors for the subsequent time levels indicate that 
the flow within the bubble is unsteady. In addition to rotating about 
the x-axis, the individual cells within the bubble tend to merge and 
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(b) t = 85.27 
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(e) t = 91.82 


Fig. 5.17 Projected velocity vectors over the interior of the 
a,b,c,d,e. breakdown region at different time levels. 


separate, and change in strength and location. The location at which 
fluid enters the bubble appears to have shifted towards the back for. the 
time levels t=89.63 and t=91.82. In addition, at these time levels, the 
forward recirculation region has lost considerable coherence. 

Particle traces are displayed in Fig. 5.18. In this figure, and 
for all remaining figures, the time level represented is t=81.28. An 
examination of the different contour plots at other time levels shows 
this time level to be representative of the solution. The nine white 
traces were started at the inflow plane from a radial position withi-n 
the rotational region of the vortex. A single blue trace was started at 
the inflow plane, but from a radial position in the irrotational region 
of the flow. The red trace was started from a position within bubble in 
the forward cell. The yellow trace was started from within the "tail" 
region of the breakdown. The white traces seem to define the general 
shape of the bubble. One of these traces enters the forward cell, 
spirals about, then exits. This seems to indicate that most of the 
fluid approaching the bubble from the front is deflected around it. The 
red trace reveals that the fluid particles in the forward cell .remain in 
the forward cell until they are forced out of the bubble. The red trace 
leaves the breakdown region from the side in the forward half of the 
bubble. The yellow trace reveals that fluid enters the breakdown region 
from the aft section of the bubble. The spiral traced out by this line 
(as it enters the bubble) is in the opposite direction to the spiral 
traced out by the blue line. The yellow line spirals about in the rear 
cell before exiting the bubble from the outer edge. Note that the red 
trace never enters the aft cell and the yellow trace never enters the 
forward cell. This indicates that the exchange of fluid between the 
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Fig. 5.18 Instantaneous streamline pattern for flow with vortex 
breakdown. 


forward and aft section of the bubble is minimal. The blue trace 
spirals about the breakdown essentially unaffected. 

Figure 5.19 is a contour plot of the axial component of velocity. 
The breakdown region is located near the center of the computational 
domain. The recirculation region is defined by the outer 0.0 level 
contour line. Immediately within this region are negative valued 

contour lines. Interior to these negative valued contour lines are 
additional positive valued contour lines. Thus, along the centerline of 
the bubble there exists a region in which no flow reversal occurs. The 
contour lines intersecting the top and bottom of the domain indicate a 
decelerating external flow (which was imposed by the boundary 
conditions) . 

A contour plot of the axial component of vorticity (scaled by the 
Rossby number) is shown in Fig. 5.20. The vorticity decreases 

continuously as the stagnation point is approached. In the forward 
portion of the breakdown region, the axial component of vorticity is 
nearly zero. This indicates that the radial gradients of the swirl 

velocity are small. In contrast, the aft section of the breakdown 

region is characterized by high levels of vorticity along the axis. 
Furthermore, the radial gradients of vorticity in this region are 
high. Note that the vorticity is negative in the region of the bubble 
corresponding to the approximate location of the aft vortex ring. 

Figures 5.21 and 5.22 are contour plots of the and co z components 
vorticity respectively. The 0.0 level contours are not displayed, as 
was the case in Figs. 5.5 and 5.6. These components are due entirely to 
perturbations of the base flow since at time t=0.0 the only non-zero 
vorticity component was u> x . The asymmetric, two celled structure of the 
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Fig. 5.19 Contours of constant axial velocity. Contour levels range 
from -0.3 to 0.9 in intervals of 0.1. 



Fig. 5.20 Contours of constant axial vorticity, w . Contour levels 

- x 

range from -.0.25 to 2.5 in intervals of 0.25.. 
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Fig. 5.22 Contours of constant vorticity. Contour levels range 
from -1.50 to 1.50 in intervals of 0.25. 
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bubble is revealed in these contours. In addition, points on the axes 
of the two vortex rings can be identified. Since the vorticity of the 
two rings is of opposite sign, they rotate in opposite directions about 
their respective axes. The maximum vorticity levels in these two rings 
is approximately equal. A possible third vortex ring can be seen in the 
outer region of the forward portion of the bubble. This ring rotates in 
the same direction as the aft most vortex ring. 

Vortex lines are shown in Fig. 5.23. The radial locations of these 
lines at the inflow plane correspond to the locations of the white 
particle traces at inflow. In the approach flow the vortex lines are 
oriented in the x-direction. The transfer of vorticity from the x to 
the y and z components takes place as the first stagnation point is 

approached. Different orientations of the vorticity vector are 

observable within the bubble. Near the centerline of the bubble, the 
vorticity vector is oriented primarily in the axial direction. Near the 
outer regions, the orientation is mostly in the circumferential 
direction. Downstream of the breakdown, the vorticity vector is 

oriented in the axial direction, but with a superimposed small 

undulation. 

1 -2 

A line plot of the integral of j w as a function of axial location 
appears in Fig. 5.24. This figure shows that the enstrophy per unit 
volume of fluid remains nearly constant in the streamwise direction 
until the breakdown region is encountered. Here, large gains in 
enstrophy are realized. Downstream of the breakdown region, the 
enstrophy level returns to the levels present upstream of breakdown. 

The volume integral of the material derivative of enstrophy is 
plotted in Fig. 5.25. The distribution of enstrophy within the 
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RATE OF CHANGE OF ENSTROPHY 


AXIAL DISTANCE 


Fig. S.24 Variation of integrated enstrophy with axial location. 



Fig. 5.25 Variation of integrated material derivative of enstrophy 
with axial location. 




breakdown region is affected significantly by both vortex stretching and 
viscous action. The effect of viscosity is to offset the gains in 
enstrophy due to stretching. The enstrophy changes due to diffusion 
across the boundaries of a specific material volume are insignificant. 
The axial location at which the enstrophy is a maximum agrees with the 
location in where the rate of change of enstrophy is zero. 

Contour lines of pressure are shown in Fig. 5.26. The minimum 

pressure at any axial location occurs along the centerline of the 

/ 

vortex, with the absolute minimum occurring at inflow. A strong adverse 
gradient is encountered by the fluid as it approaches the breakdown 
region. Within the forward part of the breakdown region, the pressure 
is nearly constant. A local minimum occurs in the aft portion of the 
bubble. This point corresponds to the location of minimum axial 
velocity. The pressure distribution is asymmetric. This may correspond 
to the orientation of the aft vortex ring. 

The pressure and axial velocity along the vortex centerline are 
plotted in Fig. 5.27. A strong adverse pressure gradient extends from 
the inflow boundary to the axial location x=3.2, corresponding to the 
first stagnation point. From x=3.2 to x=7.2 a decrease in pressure is 
accompanied by an increase in the axial velocity. The axial velocity 
then decelerates rapidly even though the pressure continues to 
decrease. Oownstream of the breakdown region, the centerline axial 
velocity is accelerated toward its free stream value. Stagnation points 
are observed at axial locations x=3.G, 5.5, 8.3. and 12.3. 

Since the external flow was decelerated in the axial direction, it 
was expected that a corresponding adverse pressure gradient would 
exist. For clarity, the pressure variation along the top computational 
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VELOCITY AND PRESSURE 


Fig. 5.26 Isobar contours: contour levels of -0.6 to 0.0 in Intervals 
of 0.1. 
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Fig. 5.27 Variation of axial velocity and pressure along vortex 
centerl ine. 




boundary (y=10.0, z*5.0) is plotted as a function of axial location in 
Fig. 5.28 using a magnified scale. A rapid change in the slope is 
apparent at x=1.5. This corresponds to the approximate axial location 

in which the imposed deceleration of the free stream axial velocity 
begins. The pressure then increases linearly until the axial location 
corresponding to the onset of breakdown is reached. Downstream from 
this location the pressure decreases and then increases. This is due to 
the curvature of the streamlines (displayed in Fig. 5.29). Experimental 
wall pressure di stributions measured by Sarpkaya [25] in a tube and vane 
apparatus behaved in a similar manner. Note that the maximum values of 
dp/dx along the computational boundary are much smaller in ragnitude 
than the values of dp/dx occurring within the vortex core. 

Figure 5.29 is a contour plot of the total pressure, q /2 + p. As 
previously mentioned, where viscous forces are insignificant and the 
flow is steady, these lines can be considered as streamlines. In this 
sense, the aft recirculation cell is clearly visible. 

The volume integral of the material derivative of energy is plotted 
in Fig. 5.30. The interpreta tion of the various lines follows from the 
description of 5.13. Note that the rate of change of energy due to 
pressure work follows the total rate of change of energy almost 
exactly. These lines appear on top of one another in Fig. 5.30. As is 
required, the rate of change of internal energy is positive throughout 
the flow. The magnitude of this change is small when compared to the 
rate of change of kinetic energy throughout most of the region. The 
plot shows that the rate of change of kinetic energy is due primarily to 
differential pressure forces acting on a material volume. The viscous 
forces responsible for changing the kinetic energy are small. 
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Fig. 5.29 Contours of total pressure: contour levels of -0.4 to 0.6 in 
Intervals of 0.1. 
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The integral of the volume averaged kinetic energy has been 
computed and plotted as a function of axial location in Fig. 5.31. The 
kinetic energy is observed to decrease by approximately 50% between the 
inflow plane and the axial location, x=8.0. Downstream of this point, 
the kinetic energy increases slightly and then remains constant. 

Axial, transverse (radial), and spanwise (swirl) velocity profiles 
at four different axial locations are plotted in Fig. 5.32 a-d for time 
level t=81.28. The axial location of the profiles in Fig. 5.32a is 
upstream of breakdown (x=0.41). The profiles in Figs. 5.32b and 5.32c 
represent axial locations within the breakdown region ( x=5 .46 and 
x=7.36, respectively). Profiles from. downstream of the breakdown region 
are plotted in Fig. 5.32d ( x= 13 .92 ) . For reference purposes, the first 
stagnation point is located at x=3.0. Upstream of breakdown, the 

spanwise (swirl) velocity profile approximates the two-dimensional 
Burgers vortex. The radial velocity is small and the axial velocity 
profile nearly uniform. Within the breakdown region, as revealed in 

Figs. 5.32b and 5.32c, the flow is no longer symmetric. Transverse 
(radial) velocities are significant and the Spanwise (swirl) velocity 
profiles no longer approximate Burgers vortex. At both axial locations 
within the breakdown, the axial velocity profiles have a local maximum 
and two local minima. The profiles in Fig. 5.32d, downstream of the 
breakdown region, reveal a large axial velocity deficit near the vortex 
centerline. The spanwise (swirl) velocity profiles reveal a solid body- 
like rotation near the centerline. At the edge of the core, the maximum 
spanwise (swirl) velocities are considerably less than values occurring 
upstream of the breakdown region. 
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Fig. 5.30 Contribution of various energy transport terms to the 
integrated rate of change of energy. 
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Fig. 5.31 Variation of Integrated kinetic energy with axial location 
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Fig. 5.32 Radial variation of axial, swirl and radial velocity 
a,b,c,d. components at different axial locations. 
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A comparison of this computed breakdown structure (at t=81.28) with 
mean streamline and axial velocity profiles constructed by Faler and 
Leibovich [27] using experimentally measured data is shown in Fig. 
5.33. The profiles given by Faler and Leibovich are for the upper half 
of the bubble only. Although the Reynolds number for the experimental 
case was of order 1000, the overall qualitative agreement regarding the 
structure of the bubble is excellent. 

Other Cases 

Several additional test cases were computed to determine the effect 
of a jet-like axial velocity profile on vortex breakdown. The results, 
although not available in the form of line and contour plots, will be 
discussed for the purposes of comparison. 

Several additional cases were computed in which the vortex was 
imbedded in a uniform external flow. For two of these cases the Rossby 
number at inflow was 0.80 and the Reynolds number, 225. The axial 
velocity profile at inflow was given by Eq. (4.6). For the two test 
cases, 6, was chosen as 1.0 and 0.5, respectively. 

The purpose of these two computations was to determine the effect 
of viscous diffusion on the flow, as measured by the local Rossby 
number. In both cases, the local Rossby number increased in the 

streamwise direction as the flow evolved. This was because the rate of 
decrease of the solid body rotation of the vortex (near the centerline) 
was greater than the rate of decrease of the axial velocity component 
(at the radius of maximum swirl velocity). The radius, r , increased 
slightly in the streamwise direction. Breakdown did not occur, since at 
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Fig. 5.33 Comparison of the numerical solution (a) at ^81.28 with 
a,b,c. experimentally determined mean streamline pattern (b) and 
mean axial velocity profiles (c) as measured by Faler and 
Leibovich (1978). 



inflow the Rossby number was above critical value and, as the flow 
evolved, it remained supercri tical everywhere. 

Two additional test cases were computed with the parameter, a, 
being chosen as 1.0 and 0.5, respectively. The Reynolds number was 225 
and the Rossby number at inflow, 0.625. The purpose of the tests was to 
determine if the jet-like axial velocity profile would delay the onset 
of vortex breakdown. The computations were not carried out to 
completion because it became apparent that breakdown would occur near 
the inflow boundary. This was analagous to the results that were 
obtained for the case Ro=0.625, Re=225 and 6 = 0.0, discussed in detail 
previously. 

Tests to ascertain the effect of grid size and domain length on the 
time evolution of the solution were performed. In all cases the Rossby 
number and Reynolds number were maintained at values of 0.5 and 225, 
respectively. The vortex was imbedded in a uniform free stream. The 
inflow axial velocity profile was uniform. The solution was computed 
over three different domain lengths; 16.0, 24.0, and 32.0 

(x-direction) . The stretching coefficient, a, -present as a parameter in 
Eq. (4.16), was equal to 1.15 for the cases h x =16.0 and h x =24.0. For 
the case h x =32 it was set equal to 1.25. This resulted in minimum cell 
sizes of 0.13, 0.195 and 0.46 for the domain lengths 16.0, 24.0, and 
32.0, respectively. The stretching parameter e, (present in Eqs. (4.17) 
and (4.18) for the y and z direction grids) was set equal to 4.5. The 
domain in the y and z direction was given by 0<y<10 and 0<z<10. In 
each case, the minimum axial velocity was plotted as a function of 
time. x results were plotted in Fig. 5.34. The results showed that 
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DOMAIN LENG 



minimum axial velocity variation with time. 



the effect of increased domain length and cell size had a minimal effect 
on the time evolution of the solution. 

It has been suggested that vortex breakdown is a consequence of 
inertial wave disturbance on the swirling flow. The azimuthal modes 
(excluding n*0) of these disturbances correspond to the asymmetric 
internal structure of the bubble. Any flow variable can be described 
by Fourier series in the azimuthal (9) direction in the form: 

l C (x,r,t) e 

- fl 

n 

It is of interest to examine the Fourier coefficients of the radial 
velocity component in front of and within the breakdown region. In 
order to compute the coefficients, the original data were interpolated 
to produce data at specified cylindrical (r,9,x) locations. The Fourier 
coefficients (C n ) were then easily computed using Fast Fourier Transform 
techniques. For these calculations, d 9 was taken as ten degrees (n/I8) 
which resulted in 36 special locations along a given radius. The 
coefficients were tabulated in magnitude, phase angle form for the modes 
n=0,l,2,3. These results, for combinations of Tx.r.t), appear in Tables 
5.1 and 5.2. The magnitudes of the coefficients for all higher mode 
numbers were insignificant, thus they were not tabulated. 

A brief investigation into the upstream propagation of wavelike 
disturbances was performed. The purpose of the tests was to determine 
if appropriately defined wavelike disturbances introduced at the outflow 
boundary would propagate upstream. Disturbances of the form 

. r 2 

5v/dy/<bw/5z = e A si n ( 2 it f * t) 


(5.3) 


Table 5.1 


Complex Fourier Coefficients for the Radial Velocity Component 


(x,r) rR3 R5I n^2 n^I 


(0.56,0.5) 

(0.56,1.0) 

(0.56,1.5) 

(0.56,2.0) 

(0.56.2.5) 

(4.15,0.5) 

(4.15.1.0) 

(4.15.1.5) 

(4.15.2.0) 

(4.15.2.5) 
(6.95,0.5) 

(6.95.1.0) 

(6.95.1.5) 

(6.95.2.0) 

(6.95.2.5) 
(8.62,0.5) 

(8.62.1.0) 

(8.62.1.5) 

(8.62.2.0) 

(8.62.2.5) 
(12.40,0.5) 

(12.40.1.0) 

(12.40.1.5) 

(12.40.2.0) 

(12.40.2.5) 


(0.03,0.0) 

(0.04,0.0) 

(0.04,0.0) 

(0.03,0.0) 

(0.03,0.0) 

(0.00,3.14) 

(0.08,0.0) 

(0.23,0.0) 

(0.24,0.0) 

(0.18,0.0) 

( 0 . 00 , 0 . 0 ) 

( 0 . 10 , 0 . 0 ) 

(0.17,0.0) 

(0.16,0.0) 

(0.16,0.0) 

(0.05,0.0) 

(0.03,3.14) 

(0.10,3.14) 

(0.14,3.14) 

(0.12,3.14) 

(0.01,3.14) 

(0.01,3.14) 

(0.02,3.14) 

(0.02,3.14) 

(0.01,3.14) 


(0.01,3.11) 
(0.01,3.09) 
(0.01,3.14) 
(0.00,3.04) 
(0.00,2.95) 
(0.00,2.94) 
(0.00,-0.78) 
(0.01,-2.29) 
(0.01,-2.77) 
(0.01,-2.76) 
(0.02,-2.92) 
(0.04,2.51) 
(0.03,1.25) 
(0.02,1.13) 
(0.01,1.35) 
(0.05,-0.46) 
(0.02,-1.17) 
(0.02,-2.06) 
(0.03,-2.27) 
(0.03,-2.06) 
(0.07 ,-1.45) 
(0.04,-1.56) 
(0.03,-1.67) 
(0.02,-1.71) 
(0.02,-1.79) 


(0.00,-1.82) 
(0.00,-1.96) 
(0.00,-1.94) 
(0.00,-1.47) 
(0.00,-1.03) 
(0.00,-0.98) 
(0.00,-0.43) 
(0.01,3.05) 
(0.00,-0.91) 
( 0 . 00 , 0 . 66 ) 
(0.00.-0.67) 
(0.02,-2.40) 
(0.02,-2.83) 
(0.01,-2.60) 
(0.00,2.52) 
(Q. 00,-0.75) 
(0.00,0.07) 
(0.00,2.50) 
(0.00,3.12) 
( 0 . 00 , 0 . 01 ) 
( 0 . 00 , 2 . 10 ) 
(0.00,1.83) 
(0.00,1.57) 
(0.00,1.81) 
(0.00,1.71) 


(0.00,-3.09) 
(0.00,-0.41) 
( 0 . 00 ,- 0 . 68 ) 
(0.00,-0.83) 
(0.00,-0.39) 
( 0 . 00 , 0 . 88 ) 
( 0 . 00 ,- 2 . 22 ) 
(0.00,-0.90) 
(0.00,-1.47) 
(0.01,-1.94) 
(0.00,2.78) 
(0.01,-1.19) 
(0.0, -1.87) 
(0.00,-0.92) 
(0.00,-0.72) 
(0.00,-1.26) 
(0.00,2.74) 
(0.00,-2.38) 
(0.00,-2.18) 
(0.00,-3.00) 
(0.00,-0.57) 
( 0 . 00 ,- 0 . 02 ) 
(0.00,1.77) 
(0.00,2.45) 
(0.00,-2.17) 
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Table 5.2 


Complex Fourier Coefficients for the Radial Velocity Component 


(x,r) n=U n^I n^T 


(0.56,0.5) 

(0.56,1.0) 

(0.56,1.5) 

(0.56,3.0) 

(0.56,2.5) 

(4.15,0.5) 

(4.15.1.0) 

(4.15.1.5) 

(4.15.2.0) 

(4.15.2.5) 
(6.95.0.5) 

(6.95.1.0) 

(6.95.1.5) 

(6.95.2.0) 

(6.95.2.5) 
(8.62,0.5) 

(8.62.1.0) 

(8.62.1.5) 

(8.62,2.0) 

(8.62.2.5) 
(12.40,0.5) 

(12.40.1.0) 

(12.40.1.5) 

(12.40.2.0) 

(12.40.2.5) 


(0.04,0.0) 

(0.05,0.0) 

(0.05,0.0) 

(0.04,0.0) 

(0.04,0.0) 

( 0 . 02 , 0 . 0 ) 

( 0 . 02 , 0 . 0 ) 

(0.19,0.0) 

(0.27,0.0) 

(0.23,0.0) 

( 0 . 01 , 0 . 0 ) 

( 0 . 01 . 0 . 0 ) 

( 0 . 00 , 0 . 0 ) 

(0.03.0.0) 

(0.13,0.0) 

(0.02,3.14) 

(0.07,0.0) 

(0.07,0.0) 

( 0 . 02 , 0 . 0 )- 

(0.01,3.14) 

(0.02,3.14) 

(0.06,3.14) 

(0.10,3.14) 

(0.11,3.14) 

(0.09,3.14) 


(0.01,3.08) 
(0.01,3.09) 
(0.00,3.04) 
(0.00,2.93) 
(0.00,2.82) 
(0.02,-2.27) 
(0.01,-2.60) 
( 0 . 01 ,- 2 . 11 ) 
(0.02,-2.29) 
(0.01,-2.28) 
(0.05,2.78) 
(0.05,-2.58) 
(0.05,-1.32) 
(0.05,-0.57 ) 
(0.01,0.61) 
(0.07,0.52) 
(0.08,1.20) 
(0.05,1.46) 
(0.04,2.86) 
(0.06,-3.00) 
(0.10,-1.76) 
(0.05,-1.97) 
(0.03,-2.35) 
(0.02,-2.57) 
(0.01,-2.42) 


(0.00,2.29) 

(0.00,2.42) 

(0.00,2.57) 

(0.00.2.79) 

(0.00,3.09) 

(0.00,-2.50) 

(0.00,-1.94) 

(0.00,-1.34) 

(0.00.-0.33) 

(0.00,0.39) 

(0.01,-1.50) 

(0.01,-0.54) 

(0.01,0.85) 

(0.02,1.77) 

(0.01,2.91) 

(0-06,0.05) 

(0.01,-0.67) 

(0.01,-2.43) 

( 0 . 01 ,- 0 . 00 ) 

(0.03,-0.13) 

(0.00,1.96) 

(0.00,1.87) 

( 0 . 00 , 2 . 86 ) 

(0.00,0.95) 

(0.00,0.48) 


(0.00,-2.39) 

(0.00,-0.45) 

(0.00,-0.79) 

(0.00,-0.92) 

(0.00,-0.92) 

(0.00,1.55) 

(0.00,0.71) 

(0.00,0.87) 

(0.00,0.29) 

(0.00,-2.54) 

(0.01,2.26) 

(0.01,-2.90) 

(0.00,-1.97) 

(0.01,-2.56) 

(0.00,-1.27) 

(0.02,1.74) 

(0.01,1.74) 

(0.01,2.60) 

(0.01,-2.85) 

(0.01,3.03) 

(0.00,1.26) 

(0.00,1.64) 

(0.00,1.89) 

(0.00,2.74) 

(0.00,-2.70) 
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were introduced where the amplitude A=0.03 and the frequency f=4.0. 
Three test cases were investigated. In two cases the vortex was 
imbedded in a uniform free stream. In one of these cases the Rcssby 
number was 0.5 and in the other case the Rossby number was 0.8. In the 
third case the- Rossby number at inflow was 0.8 and the vortex was 
imbedded in a decelerating flow. In each case the solution obtained 
with forcing at the outflow boundary was subtracted from a base flow 
solution. This base flow was the solution that resulted from identical 
inflow and initial conditions, but without forcing at outflow. For the 
specific forcing frequency and forcing amplitude combination applied, no 
disturbances were seen to propagate upstream (for time t<8). A 
discussion of the motivation for these tests follows. 

The theories of Squire [7j and Benjamin [8] pertain to stationary 
disturbances present on columnar vortices. Randall and Leibovich [18] 
have shown that for upstream propagation to be possible, the base 
vortical flow must change in the axial direction. This would be the 
case if an adverse pressure gradient were encountered. This pressure 
gradient could be self-induced (through vi scous- di ffusion) or externally 
imposed. Thus the motivation for the three previously mentioned 
tests. That the disturbances did not propagate upstream in the cases 
where the vortex was imbedded in a uniform free stream was expected. 
These solutions were allowed to develop for 8 time units, thus the 
vortices remained nearly cylindrical (no variations were present in the 
axial direction). In the third case it was expected that some type of 
disturbance would have propagated upstream to the critical station. 
Possibly the forcing frequency /ampl i tude combination that was imposed 
resulted in waves that were damped immediately. 


110 


CONCLUSIONS 


Numerical solutions of the fully three-dimensional Navier-Stokes 
equations were obtained for vortex breakdown. The numerical algorithm 
was an implementation of the veloci ty-vortici ty formulation developed by 
Gatski, Grosch and Rose [39]. In this formulation, both the velocity 

and vorticity vectors are second order accurate in space and time. The 
solutions were presented for unconfined vortices of the leading edge 

type and were parameterized by the Rossby number and the Reynolds 

number. Breakdown was predictable using the local Rossby number as the 
critical parameter. 

The present analysis supports Squire's [7] earlier conjecture that 
the physical mechanism responsible for breakdown is the growth of 
wavelike disturbances along the vortex core. The Coriolis acceleration 
produces a restoring force which is responsible for the creation of the 
disturbance waves. This effect was described in terms of the Rossby 
number by examining the theoretical, computational and experimental 
resul ts available in the literature. The resulting correlations showed 
that when the Rossby number of the base vortical flow decreased to 

approximately 0.6, breakdown occurred. 

The Rossby number criterion may find practical applications in the 
field of aeronautics. For instance, it would be advantageous if vortex 
breakdown could be induced in the case of trailing wing tip vortices 
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generated by commercial aircraft. Retarding vortex breakdown is 
desirable in the case of leading edge vortices generated by delta wing 
type military aircraft. In applications such as swirl combustors, the 
internal structure of the breakdown region is of importance. The size, 
shape and stability of the recirculation zone are critical to flame 
stability and performance. Optimal rates of entrainment and mixing 
could be predicted through a numerical simulation of the process. 

Detailed results of vortex breakdown were obtained for two 
numerical simulations. In the first case, (Ro=0.625, Re=225) the vortex 
was imbedded in a uniform free stream. In the second case, (Ro=0.8, 
Re=225) the vortex was imbedded in a decelerating free stream. The 
internal structure of the resulting solutions differed dramatically. 

The structure resulting from the first case consisted of a single, 
steady, nearly symmetric toroidal recirculation zone. The length of the 
bubble was defined- by the two existing stagnation points. The structure 
in the second, decelerating case contained multiple internal cells, or 
"ring" vortices. In addition, the flow within the bubble was unsteady, 
asymmetric and dominated in the aft portion by 'large amplitude velocity 
fluctuations. 

The discrepancy between the results of the first case and the 
experimental results may be attributed to the proximity of the numerical 
breakdown to the inflow boundary. In that case, the fixed boundary 
conditions at inflow act as an artifical critical barrier to the 
propagation of waves. At the inflow boundary, the amplitude of a 
standing wave grows to the extent that stagnation points and reversed 
flow appear along the axis. In the decelerating case, the critical 
condition arises from the Rossby number effect in a manner analagous to 
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experimental findings. The solution is not contaminated by the inflow 
boundary layer or by an artifical critical condition. In addition, the 
unsteady features of this case may be due to disturbances originating 
downstream since the flow was imbedded in the decelerating stream. The 
bubble in the first case was imbedded in a uniform free stream, thus it 
is expected that downstream disturbances should be damped immediately. 
This would account for the nearly steady features of that flow. This 
could not be confirmed by the tests which employed forcing at the 
outflow boundaries. 

For the case Ro=0.8, Re=225 the radial velocity component was 
expanded in a complex Fourier series in the circumferential direction. 
An examination of the resulting Fourier coefficients c n (x,r,t) has shown 
that the n=0 mode is the dominant mode at a location slightly upstream 
of the first stagnation point for all time levels. At time level 
t*8 1.28, the magnitude of the n=l mode in certain areas within the 
breakdown region is significant. At the later time level t=91.82, the 
magnitude of the n=l modes within the breakdown region has increased. 
This indicates that the asymmetric structural- features of the bubble, 
due primarily to n=l mode disturbances, have been enhanced with time. 

In conclusion, theoretical results, interpreted in terms of a 
Rossby number, give conditions favorable for the occurrence of 

breakdown. By considering the flow to be three-dimensional, unsteady 
and asymmetric, the author believes the first correct numerical 

representation of the internal structure of the ‘‘bubble" type breakdown 
was computed. 
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APPENDIX A 


EXPONENTIAL TRANSFORMATION 


Gatski, Grosch and Rose [39] employ an exponential transformation 


of the form 


i [B (t-t )] _i 
u = e n n C 


where [B n ] is a constant matrix, to eliminate the explicit appearance of 
the vortex stretching term from the vorticity transport equation. The 
procedure is outlined below. 

The vorticity transport equation is written in the form 

lr* u j-V“ Jb V' j -« (a2) 

. B (t-t ) . 

let u = B n , a constant at time level n. Introduce w 7 = e n n C 1 

into the above equation. This results in the following transformed 


equa tion. 


[8 (t-t)] . [B ( t-t )] ft _i [B { t-t ) ] 

B n e n " c’ + e " n -$■ * e n n Uj c’.j 


i CSn^-r)] , 
- 3 n e " n C J = ve " " c' 


It is apparent that the vortex .stretching term is eliminated 
through the time derivative. The resulting equation appears below. 


§§- + u. C 7 , • = vC 7 ,.. 

at J J JJ 
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APPENDIX B 


REDUCTION OF X-SWEEP VORTICITY EQUATIONS TO 
TRIDIAGONAL FORM 

The reduction of the x- sweep vorticity equations to tri diagonal 
form is shown below. For the sake of clarity only appropriate indices 
are expressed. The derivation of the y and z sweeps are analogous, with 
only the final results being expressed. 

The 3-D advection diffusion equation for the transformed vorticity 
and the necessary averaging conditions appears below as Eqs. (Bl) to 
( B14 ) . Note that the (u,v,w) velocity components represent averages at 
cell centers. 


+ u6^C + v ^yC + w§ z C 55 v(5 x $ + (Bl) 

<“x - Wx> ♦ * 6 x" C ' (B2) 

S ' VyV * * V (B3) 

! “z ’ h z , x 6 z ) 1 (84) 

“t 5 * p x’ c ' ” 2 x p x 6 x* (85) 

■ 4 y C - h 2 y P y 6 y * (36) 
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where 



^ t C « n 2 C - 

" 2 p z 6 z 5 


Ax 

(B8 ) 

n a 

q<V 

= 2 

p y 

9 

y 

u 

1$ 

(B9) 

n s 

q( 9 z ) 



p z 


AZ 

2 

( BIO) 

and. 

z 

At 



T 3 “I 

q( 9 ) 




X 

" 6 

( Bll) 




( B7 ) 


(B12) 

(B13) 
(B14) . 


The first step is to eliminate the n+ V 2 time level dependence from 
equation (Bl) as follows. Combine the identity 

_ - n- V 2 


x6 t C = n t C - C 


(B15 ) 


with Eqs. (B5) 
equations. 


through (7) to obtain the following three vector 

2 - n- V? 

x6C = u C - h p 6 * - C 2 (B16) 

t X XXX 

— — 2 ** ** n- V? 

■c6 C 8 (i C * n p 5 <I»~C 2 ( B17 ) 

t y y y y 

2 - n- V? 

tfi C * n C - h p 6 5 * C 2 (B18) 

u u Z Z Z 


When Eqs. (B16) to (B18) are each combined with Eq. (Bl), three 

systems of equations result, each of which can be solved simultaneously 

for a component of the viscous terms 6 $, 6 4 * and 6 

x y z 

The systems appear below as Eq. (B20). 


K 


x 


1 


1 


vx 


1 


K 


y 


1 


1 


1 


K 

z 


6 $ 
x 


5 4 = T 

y 



1 1 1 
1 1 1 
1 1 1 


u6 x c 

v6 y C 
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where 



0 0 
1 0 

0 1 


\i c - 
x 



^ ’ 


l V 2 
l n_1/ 2 
C "* V 2 


h 2 p 

K 1 

X VT 


(B19) 


(B20) 


h 2 P 

k = + i 

y vt 


(B21) 


t, 2 P 

K = — i <• 1 

2 VT 


(B22) 


Take the inverse of the matrix on the left side of Eq. (B19). This 
is given as 


where 


11 

a 12 

“l3 

21 

a 22 

“23 

31 

“32 

“33 


(B23 ) 


M 2 + K x K y K z - (K x + K y + K z > (B24) 


11 * 

(K K - 1) 
y 2 

( B25) 

12 = 

- (k z - 1) 

(B26) 

13 = 

(1 - K y ) 

(B27 ) 

21 = 

-(k z - 1) 

( B28 ) 

22 * 

(k x k 2 - 1) 

( B29 ) 


19 ? 



“23 ’ "< K x ' 11 


(830) 


“31 * (1 ' V 


“32 ■ * (K x ‘ » 


“33 * (K x K y - 11 


Multiply matrix Eq. (B19) by matrix (B23). This results 
following expressions for the viscous terms. 


1 f , - - n- V? . 

= - — a., u C • C c ) + 


54 ,= - — { a ( \i C - C ) + a ( p. C - C ) 
x 11 x 12 y 


* “13 ‘“z'* ' "" 1/2)1 


+ ^ (a ll + a 12 + a 13 )(u6 x C + v 6 y C + w 6 z C) 


V = f “2l‘ W x' ‘ ' ^ 1 * “22V ' E V2 1 


* “23 ( V“ " c " 2)1 


+ <“21 * “22 + “23 Hu6 x C + , 6 y C + " 6 z° 


1 t s (“-,,<1* c - C " ) + a (n C * C n ^ ) 

z Bvt 1 31 x 32 y 


+ “ 33 U z i ’ 5 "' 1/2)1 


"vS (a 31 + a 32 + a 33 )(u6 x’ C + v6 y~ C + w6 z^ 


(B31) 
( B32) 
(B33) 
i n the 


(834) 


(B35) 


(B36) 
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Only Eq. (B34) is needed for the derivation’ of the x-sweep tridiagional 
system. Equations (B35] and (B36) are needed for the y-sweep and z- 
sweep derivations respectively. Performing the difference operation 
indicated on the left side of Eq. (B34) results in the following 
expression 


U. i.) = -r— - (a ( m. C - C ^ ) 

1+V2 l - 72 pv\^ 11 * 


+ a i2 (4 y^ ’ ^ 0 + “l3^z^ " ^ 0 ^ ^ 



Ax , . , 

+ — (a + a + a ) ( u 6 C + v6 C + w6 C) 
0v 11 12 13 x y Z 

{ B37 ) 

where 

X 

II 

• 

(B38 ) 

Let 

Y x ± 3 V 2 (l ± q( 9 X ) ) 

(B39) 


The appropriate combination of Eq. (B2) and Eq. (B39) results in 
the following expression for y 

\ v * 6 x" c • \ *i- v 2 • • (B40) 


Multiply each side of Eq. ( B37 ) by y* and substitute the result into Eq. 
(B40). The resulting expression for + y is given as. 


2y 


^i* 1/2 ~ v 11 x 


(a (u C * C " ) +* a ( u C ” C 0 ^ ) 

12 y 


C 

13 z 


- C 1/2 ) 1 
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( B41 ) 


+ 

Uy Ax 

+ [ ~V- ( “ll * °12 * ‘l3> + ^ V 


x Ax 


ve (a ii* e ‘i2 t “i3 lc,5 y ,: + “ 6 2 ,:] 


The expression for y , obtained in a similar manner, appears 


below. 




* «13 ' C n " V2 )l 


-Uy 


x Ax 


.p- ( “ll * “12 + “l3> + '] *,< 


Tx + Ax 


v? ( “ll * “12 + “l3 )Cv V * w6 x=l 


(B42) 


Index Eq. (B42) by i+1 and equate fluxes of <p across a cell face to 
obtain the following equality. 


- c « u (1 X C 


1+1 ■ c a ll u x c 


+ [- UR- (« u + u 12 « 13 ) + l]6 x C 


i+1 


- [ uR + a i2 a 13 ^ + ^ C 


■ c « 12 v-'"' V2 » 


i+1 
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* c* « (» c - C "* 1/2 

12 y 


4 C ' “13 ' « n ' 1/2 » 


1+1 






+ - n- 1 h 

1+1 * C “ll C 


+ R (<x n + a 12 + « 13 ) [v6 y l + w6 z l] 


1+1 


+ R (<x u + + a 13 )[v6 y C + *5^] 


where 


.* V 

• Sv\ 


and 


0 ± X AX 
R " 0v 


Performing the discretizations, rearranging and redefining y : 

T x* 1 1 4 o' s x 1 

results in the tridiagonal form of the equations for the x sweep 
below. 


-Nil K (-f'-TT' ' s x ] |--- c - 


a a uX 
11 1 x. 


x | i +1 i+ 3/2 


+ Nil 


a,. a u\ 




( B43) 

(B44) 

( B45 ) 
as 

( B46) 
given 



,i — ^*1 1 a 1 u\ 


. .-1 r * , V . V x 


X ?[*x «-r*-V-- s ,]|i?i-W 


Zy a. 
x 1 


\ 0 

X 


1 / 2y a. 
C V 2 - _i 1 


1 + 1 


1 + 1 


\ 0 
X 


- n- V 2 


C. 

i 


* r * ,“12 . VY 

x ( t ~r L) 


i+l C i+l,j+V 2 


♦ :*_££. vv 

s 8 8 


i+l C 1+1. j- V 2 


Y v a.- a.vX 

+ 4- (-P + _J Z) 

\ * p 


1 C i.j+V? 


+ 1* (li! 

X x P 


a, v\ 

1 ,y) 
s 




+ I*. fl3 . VS, 

\ 8 8 


i+l C i+l,k+l/ 2 


* I*_ ( fl3 VY 

\ 8 8 


i+l C i+l,k-V 2 


r x “u “i"S 

+ _i_ (_ii + i) 

X x 8 8 


f Yk+V 2 
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+ I* ( !l 3 


. a i w v 

ft 3 1 


i “*.k-V 2 


where 


“l = a ll + a 12 + “13 


and 


2v\ 


= 


x Ax 


The tri diagonal equation for the y sweep takes the form 


a oo auvX. 

j+1 V ( ~F + " V|j+i 


- 


j+ 3/2 


* t T y * S y]|j +1 


— 1 - 4. ^99 (Z*)V\ 


1 . ^99 Cl<jV \ 

’ x j Y y ( 1 r^ } “ S yJ| j \j- 


V? 


• 2y y “2 


V 


1/ a Y a >, 
- n- V 2 _ y 2 


c . 


i+i j+1 *• e 

J+l y 


- n- y 2 

j 


+i-<‘ 2 i+i 


U V. 


T* ( T 

y 


j+i ; i + V 2 , j+1 


, T y ,“21 “2 U \ 

X* ( "5 5 -' 


j+i V 2 ,j+i 


( B47 ) 
(B48) 
(B49) 
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Y y *21 OoUX 

* 4- Hr + -V' 

y H p 


i c i ♦ V 2 ,j 


i (“21 . “2 U \, 

\ 1 9 ~1~' 


i C i'-V 2 .j 


i. (“23 + a 2 w S 


X 8 

y p 


3 


j,( 'j + l.k+V 2 


* lK_ (^23 - V^z, 

\ 9 9 ' 


j+1 C J*l.k-V 2 


r v a 22 <z ? w '^7 

+ -f 


j c j, 1 <»v 2 


,IL. ,“23 + “2“\ 

* -i- ( -r + -5-) 


j c j, k- V 2 


Where a 2 = + a 22 ♦ ^3 


and 


2vX^ 

s y " “^7 


The tridiagonal equation for the z-sweep takes the form 


-1 - a V3 CC-WX 

X k+1 t y z { ~ + F“ } “ s z^|k+l ^k+ 3/2 


, _i a,, a-wX 

+ { " H+i t y z { “r " + s z^|k+i 


(B50) 

(B51) 

(B52) 
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<-P* 3 z 


z ( 8 + f 1 + S zl|kl ^k +IA 


i ^*3^ CX-WX. 

"1 r * / 33 3 z x 


X l^ T z ( 9 p 1 s z J|k 5 k _ l/ 2 


2y z *3 


V 


- n- V 2 _ 2y z “3 


k+1 


k+1 


V 


r n " :/ 2 

k 


* 2 0 • 0 ' 


k+1 


*1+ V2 ,k+l 


+ li. ,^31 . VS. 

\ 'a a I 




3 


k+1 ^2 »^ + l 


y „ a a u\ 

+ — (— + J. . X ) 
* 2 0 0 


’i+ V2 ,k 


+ Y z .“31 .V\c. 

\ l “0 F~ } 


k C i" V 2 ,k 


Is f “32 . V»» } 

V 0 0 J 


k+1 C J +1 /2,k+l 




k+1. 



V 2 ,k+l 


i ^?n 


T 


+ 



22 

3 


vX 


c j 


*V 2 .k 


+ 





( B53) 


where 


= a. 


31 


+ a 


32 + a 33 


( B54) 


and 


s z * 


2v\ 

2 

AZ 


(B55) 
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APPENDIX C 


DERIVATION OF PRESSURE EQUATION 


The three scalar components of the Navier-Stokes equations are 


given by 


A au ^ au ^ au -1 ap 3. 

—T + u -7— + V — - + W - = r— + vv U 

at ax ay az p ax 


av , av , av , av -1 ap A _2 

at ax ay az p ay 


aw aw aw aw -i ap _2 

-r— + u — + v-r- + w— = -~rr- + vv w (C3) 

at ax ay az p az 

Take the partial derivative of Eq. (Cl) with respect to x, Eq. (C2) 
with respect to y and Eq. (C3) with respect to z and find 


(£U) + (®it) + u a«i + v -LiL 

at 'ax' 'ax J ax ay - axay 


aw au a u _ a p a , _2 > 
+ " aJH + * aJaz T77 + v aT (v u) 

OX 


2 2 2 

a , av» au av a v , av, , .. a v L aw av 

at ( a?’ a? aJ * “ aSa? + ‘a? 1 + v ^ * a7 a? 


, , a 2 v l a 2 p ^ a ,„2 , 

ayaz p gy2 ay 


d_ ,dw, + 3 2 * + aw a 2 w . dw. 2 

at az dz dx dxdz az ay ayaz az^ 


2 ? 

, a w l a“p a , 2 > 

77*'^ + u « 17 


Add Eqs. (C4), (C5) and (C6) together to find. 


±- (*H + *1 + aw) ( au . 2 ,av 2 , aw 2 av au 

at l ax ay az ; W + ( a? ( a? + 2 ( a? ay 


,2 .2 


aw au aw av. ,au. a v aw . 

ax az a y az } u ( 7T axay aTal* 


+ V {iz + £lL ♦ ®?!L + w ( ifiL + SJL + s 2 w. 

axay ayaz w l axa y aPal + ~ 1 ) 


= -i7 2 p + vC|-(7 2 u) ^I-IT 2 .) *|-(7 2 „)] 


Imposing the continuity condition 


au + av + aw 
ax ay az 


allows Eq. (C7) to be written as: 


2 2 2 

(1H) + (SV) + + 2 f— + *S du + Sw av. _ -l ? 

W V ( az J + 2( a7 ay + lx al + ay az } T p< 


Expand the first three terms on the left side of Eq. (C9) as 


(2£) + (±1) 2 + = ( au + av aw 2 

ox 'ay az ax ay az 


(5u av + au aw + av au + av aw + aw au aw av, 

ax ay ax az ay ax ay az az ax az ay 
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Imposing the continuity equation on Eq. (CIO) and substituting the 
results into Eq. (C9) results in the following elliptic equation for the 
pressure field. 

_2 ? r£v au . au aw au av au aw av dw . 
p '^ PL ax ay ax az ay az ax ay ” ax az ay az J lcllJ 
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